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Preface 

This  report  gives  somewhat  abbreviated  versions  of  talks 
in  a  seminar  on  geophysics  that  was  held  once  a  week  during  the 
academic  year  1970-71  at  the  Courant  Institute  of  Mathematical 
Sciences  of  New  York  University. 

The  talks  were  presented  for  the  most  part  by  applied 
mathematicians  rather  than  by  professionals  in  the  field  of 
geophysics,  but  a  number  of  professionals  attended  the  seminar. 
It  was  hoped  that  both  groups  would  find  it  interesting  and 
rewarding  to  look  at  the  problems  together,  based  on  the  well- 
founded  experience  in  science  that  progress  has  often  resulted 
by  bringing  disciplines  together  that  at  first  sight  might  not 
seem  to  have  much  in  common.   It  was  also  thought  desirable  to 
examine  a  rather  wide  variety  of  physical  problems,  which  in 
turn  required  the  use  of  a  variety  of  mathematical  ideas  and 
methods.   For  example,  more  emphasis  was  put  on  problems  with 
formulations  from  the  theory  of  elasticity  in  its  linear  version 
than  is  usual  in  the  literature  on  geophysics. 

Some  of  the  ideas  advanced  may  well  be  found  highly 
speculative  by  professionals,  but  it  is  hoped  that  they  may 
nevertheless  be  found  amusing  and  stimulating.   Of  course,  it 
is  also  hoped  that  applied  mathematicians  may  find  it  worth- 
while to  become  acquainted  with  the  highly  attractive  and  varied 
field  of  geophysics. 
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Lecture  1 

The  Gyroscopic  Motions  of  the  Earth 
J.  J.  Stoker 

A  variety  of  phenomena  in  geophysics,  such  as  the  diurnal 
variation  of  the  tides  in  the  ocean,  and  seasonal  changes  in  the 
atmosphere,  are  caused  by  gyroscopic  motions  of  the  earth,  which 
in  turn  are  the  result  of  gravitational  act±on  by  the  sun  and  the 
moon.   In  this  first  lecture  the  salient  features  of  the  gyro- 
scopic motions  of  the  earth  are  described. 

In  geocentric  coordinates,  the  earth  spins  with  a  period 
T   =1  day,  the  sun  moves  about  the  earth  with  a  period 
T-,  =  1  year  and  the  moon  moves  about  the  earth  with  a  period 
T      =   27.2  days.   The  earth's  spin  axis  rotates  with  a  retrograde 
precession  about  the  normal  to  the  ecliptic  (the  plane  of  the 
sun's  orbit)  with  a  period  T^  =  26,000  years,  whereas  the  moon's 

node  (the  line  in  which  the  plane  containing  the  moon's  orbit 

2 
intersects  the  ecliptic)  precesses  with  a  period  T^  =  18—  years. 

Finally,  the  earth's  axis  executes  small  oscillations  of  amplitude 

measured  in  seconds  of  arc,  called  nutations,  which  have  various 

periods  and  which  have  various  causes. 

A  description  of  these  motions  based  on  rigid  body  dynamics 
will  be  given,  but  without  detailed  analysis.   For  such  an  analysis 
the  book  of  Klein-Sommerfeld:   Die  Theorie  des  Kreisels,  should  be 
consulted. 

The  orbital  motions  of  the  sun  and  the  moon  relative  to  the 
earth  come  from  the  gravitational  forces  between  them.   The 
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gyroscopic  motions  of  the  earth  arise  from  the  gravitational 
torques  exerted  on  it  by  both  the  sun  and  the  moon.   These  torques 
owe  their  existence  to  the  departure  of  the  earth  from  sphericity. 
In  fact  the  ellipsoidal  earth  with  mass  m  has  an  axial  moment  of 
inertia  C  larger  than  its  transverse  moment  of  inertia  A  such  that 
(C-A)/C  ~  1/305.   The  earth's  radius  is  small  compared  with  the 
distances  to  the  sun  and  the  moon.   As  a  consequence  the  gravita- 
tional potential  energy  of  the  ellipsoidal  earth  may  be  replaced 
with  good  accuracy  by  that  of  a  homogeneous  sphere  plus  a  thin  uni- 
form ring  around  its  equator  in  such  a  way  that  the  total  mass,  m  , 
and  the  moments  of  inertia  of  this  combination  are  the  same  as 

those  of  the  earth.   The  ring  must  have  a  radius  r  and  a  mass  m, 

o  o 

and  the  sphere  must  have  a  mass  M,  such  that  M+m  =  m  , 

^Mr^  +mr^  =   C   and  ^Mr^  + -imr^   =  A,    i.e.      M  =  5(2A-C)/r^   and 
500  5020'  -^  ^  "    o 

m  =  2(C-A)/r  .   As  the  sun  moves  in  an  orbit  which  is  very  nearly 
circular  and  hence  with  constant  velocity,  its  gravitational 
potential  is  a  periodic  function  and  can  be  expanded  easily  in  a 
Fourier  series.   The  time-independent  secular  term  in  the  Fourier 
series  can  be  interpreted  with  good  accuracy  as  the  potential  due 
to  a  uniform  thin  ring  with  a  mass  m-,  equal  to  the  sun's  mass  and 
a  radius  r-,  equal  to  the  average  distance  between  the  sun  and  the 
earth.   Similarly,  the  moon  may  be  replaced  by  a  uniform  thin  ring 
with  a  mass  m^  equal  to  the  moon's  mass  and  a  radius  rp  equal  to 
the  average  distance  between  the  moon  and  the  earth.   These  approxi- 
mations can  afterwards  be  improved  by  considering  the  effects  of 
eccentricity  of  the  orbits  and  of  nonuniform  velocities  as  small 
perturbations . 
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The  sun's  ring  lies  in  the  ecliptic,  but  the  earth's  ring 

is  inclined  to  the  ecliptic  at  an  angle  9   approximately  equal  to 

1° 
25^  .   Thus  the  gravitational  attraction  of  the  sun's  ring  on  the 

earth's  ring  produces  a  torque  about  the  earth's  node  (the  line 
in  which  the  plane  containing  the  earth's  ring  intersects  the 
ecliptic).   A  similar  torque  is  produced  by  the  gravitational 
attraction  due  to  the  moon's  ring.   Neglecting  the  small  inclina- 
tion angle  (approximately  5°)  between  the  moon's  ring  and  the 
ecliptic,  the  resultant  torque  on  the  earth's  ring  is  found  to  be 
given  approximately  by 

^   o /m,    m  \ 

■^mr   -i  +  -1  cos  e   sin  Q  '-  ■  '      '  ■ 

where  G  is  the  gravitational  constant.   This  torque,  acting  on  the 
earth's  ring,  which  spins  vn.th  the  period  of  T  ,  produces  a  pre- 
cession of  the  earth's  axis.   The  angular  velocity  of  the  pre- 
cession is 

2 

-,      rnrT^  /   m,    m^  . 


\  r-i        r 


It  is  vjorth  mentioning  that  the  contribution  due  to  the  moon's 
attraction  is  larger  than  that  due  to  the  sun  by  a.   factor  of  2.13 
because  the  distance  to  the  moon  j. s  small  enough  to  outweigh  its 
much  smaller  mass.   (This  is  also  the  reason  why  the  moon's  effect 
is  predominant  in  creating  the  tides  in  the  oceans. )   Using  the 
relationships 
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as  given  by  Kepler's  laws  we  obtain  the  period  of  the  earth's 
precession  as 


i  =  1  C-A 
T^   2   C 


T  T^ 

o      ,      o 


cos  9 


If  the  value  of  (C-A)/C  is  taken  as  1/305,  we  get  T,  equal  to 

1" 
26,000  years  (i.e.   5I—  of  arc  per  year).   In  fact,  however,  the 

above  formula  was  really  used  to  calculate  (C-A)/C,  by  using  the 
observed  value  of  T  .   The  result  is  in  agreement  with  other 
methods.   Moreover,  the  precession  of  the  earth's  axis  is  retro- 
grade, i.e.  the  equinoxes  (the  points  where  the  earth's  node 
intersects  the  sun's  ring)  move  in  a  sense  opposite  to  that  of 
the  earth's  spin.   This  is  so  because  the  torques  tend  to  pull 
the  earth's  ring  into  the  ecliptic. 

In  the  same  manner,  due  to  its  inclination  to  the  ecliptic, 
the  rotating  moon's  ring  precesses  under  the  gravitational  torque 
produced  by  the  sun's  ring  (the  torque  produced  by  the  earth's 
ring  is  negligible).   The  calculated  period  is  18  years  under  the 

assumption  made  here,  which  is  considerably  smaller  than  the 

2 
observed  value  of  18=-  years.   The  discrepancy  arises  in  the  main 

from  the  rather  large  eccentricity  of  the  moon's  orbit. 

The  steady  precession  of  the  earth's  axis  is  subject  to 

various  small  amplitude  nutations.   The  component  due  to  the  free 
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oscillation  (as  the  initial  conditions  do  not  fit  the  exact  pre- 
cession) has  an  amplitude  of  the  order  of  tenths  of  a  second  of 
arc  and  a  period  of  504  days,  which  is  much  smaller  than  the 
observed  value  of  428  days.   The  discrepancy  is  attributed  to  the 
elasticity  of  the  earth.   It  appears  that  the  earth  should  have 
the  rigidity  of  steel  to  reconcile  the  discrepancy.   More 
important  components,  because  of  larger  amplitudes,  come  from 
the  forced  oscillations  due  to  external  periodic  torques  with 
periods  equal  to  those  of  the  precession  of  the  moon's  ring,  and 
the  fact  that  the  orbits  of  the  sun  and  the  moon  are  not  circular. 
These  terms  arise  from  the  harmonic  components  in  the  Fourier 
series  for  the  time-dependent  gravitational  potentials  of  the 
precessing  moon's  ring  and  the  orbiting  sun  and  moon. 

To  summarize,  we  combine  the  various  gyroscopic  motions 
linearly  in  the  following  expressions  (time  t  is  measured  in  years) 
for  the  inclination  angle  6    and  the  azimuthal  angle  f,    both  with 
reference  to  the  ecliptic,  which  specify  the  orientation  of  the 
earth's  axis. 

'Z'  =  50". 371^  t  -  17". 251  sin  -p-  +   0".207  sin  -J- 

T^  ill 

-  1".269  sin  ^  -  0".204  sin  -J-^  , 
^1  ^2 

0  =  23°27'32".0  +  9". 223  cos  -^  -  0".090  cos  -J^ 

+  0".551  cos  M  +  o".089  cos  ^^  . 
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The  perturbation  due  to  the  second  harmonic  of  the  precession  of 
the  moon's  ring  is  included  here.   More  correction  terms  can  be 
added  to  account  for  the  ellipticity  of  the  orbits,  motions  of 
the  perihelia, and  interaction  of  other  planets. 
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Lecture  2 


Flood  Waves  in  Rivers  and  Reservoirs 
J.  J.  Stoker 

Flow  problems  in  rivers  and  reservoirs  were  treated  in  the 
past  by  hydraulic  engineers  using  the  flood-routing  methodo'- 
Owing  to  its  oversimplification  of  the  dynamic  behavior  of  the 
flows,  the  method  was  not  always  applicable.   With  the  advent  of 
high  speed  digital  computers  it  became  feasible  to  solve  the  flow 
problems  numerically  by  finding  the  solutions  for  the  governing 
nonlinear  partial  differential  equations.   The  first  attempts 
were  made  by  Stoker,  Isaacson  and  Troesch'-    in  1953-1956  on  a 
UNIVAC  computer  (see  also  [3]).   Three  test  cases  were  made,  with 
the  object  of  comparing  the  observed  flood  stages  with  those 
computed  numerically:  ■■,         • 

(1)  Flow  in  a  long  river:   the  flood  of  19'^5  in  the  375-mile 
stretch  of  the  Ohio  River  between  Wheeling,  W.  Va.,  and  Cincinnati, 
Ohio, 

(2)  Flow  through  a  junction  of  two  rivers:   the  flood  of  19-^7 
through  the  junction  of  the  Ohio  and  Mississippi  Rivers  at  Cairo, 
111.,  about  35  miles  in  each  branch,  and 

(3)  Flow  in  a  long  reservoir:   the  flood  of  1950  through 
Kentucky  Reservoir,  about  184  miles  long,  in  the  Tennessee  River. 

In  each  case,  the  flood  stages  (the  height  of  the  water 
surface  above  sea  level)  and  the  flow  velocities  (averaged  over 
the  cross  section)  were  assijmed  known  everywhere  at  some  arbitrar- 
ily selected  initial  time.   For  subsequent  times  the  inflow  from 
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tributaries  and  the  local  run-off  in  the  main  valley  were  taken 
from  the  actual  records,  and  the  progress  of  the  flood  in  the 
main  stream  at  various  locations  along  the  river  was  then  computed 
for  future  times  of  the  order  of  weeks.   The  numerically  calcula- 
ted flood  stages  were  then  compared  with  those  actually  observed 
along  the  river.   The  results  as  shown  in  Figures  1,    2,    and  3  were 
very  satisfactory.   Thus,  it  was  demonstrated  that  numerical  flood 
predictions  have  become  possible  in  a  practical  way  (in  contrast 
with  the  earlier  failure  in  the  much  more  complicated  case  of 
numerical  weather  predictions  )  as  it  is  conceivable  that  it  might 
take  longer  to  predict  a  flood  numerically  than  the  duration  of 
the  flood  itself.   Even  on  the  UNIVAC  the  computing  time  needed 
to  compute  flows  for  periods  of  three  weeks  was  about  six  hours; 
with  the  equipment  available  now,  the  CDC  6600,  the  time  would  be 
reduced  to  seconds.   After  some  initial  hesitation,  hydraulics 
engineers  have  accepted  the  method  as  practical  and  useful,  and 
it  has  superceded  the  building  of  very  expensive  hydraulic  models 
—  which,  in  any  case,  were  not  true  models  but  rather  analogue 
computers. 

The  laws  of  conservation  of  mass  and  momentum  for  flow  in 
an  open  channel  are  formulated  in  the  following  two  partial 
differential  equations  for  H  (the  flood  stage,  e.g.  the  height 
above  sea  level)  and  u  (the  flow  velocity)  in  terms  of  the  distance 
X  along  the  river  and  the  time  t: 


* 

At  the  present  time,  numerical  prediction  of  various  quantities 

important  in  weather  forecasting  is  carried  out  successfully,  at 
least  for  times  of  the  order  of  days. 
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B(x,H)  -^  H  +  -|^[A(x,H)u]  =q(x,t)  , 

-^u  +  u-|^u  +  g-|^H  =  -F(x,H)u|u|  . 
In  these  equations  g  is  the  acceleration  of  gravity,  A  is  the 

Sa 

cross-section  area  of  the  river,  B  (=  -^)    is  the  breadth  of  the 
river  at  its  surface,  and  F  is  the  friction  coefficient.   As 
indicated  A,  B,  and  F  are  functions  of  x  and  H.,   In  any  given 
case  they  are  quantities  that  must  be  determined  empirically;  the 
first  two,  of  a  geometric  nature,  were  obtained  from  topographic 
maps  of  the  river  valley,  while  the  friction  coefficient  was 
determined  by  using  flow  data  for  past  floods  from  F  =  g|-g5^l/'^  , 
by  neglecting  the  inertia  term  in  the  equation  of  motiono   The 
source  term  q  represents  the  rate  of  the  inflow  into  the  main 
channel  per  unit  length  of  the  river,  either  from  tributaries  or 
by  the  flow  over  its  bank  from  the  main  valley;  this  quantity  was 
estimated  from  actual  records  in  our  test  cases.   The  term  Fu|u| 
represents  a  resistance  which  is  in  the  direction  opposite  to  the 

flow. 

In  addition  to  the  differential  equations,  it  is  necessary 
to  prescribe  initial  data  and  boundary  data  to  have  a  completely 
formulated  problem  with  a  definitely  determined  solution.   As 
initial  data,  the  values  of  the  stage  and  velocity  at  some  instant 
of  time  were  taken  from  flood  records.   As  the  boundary  data  for 
the  first  test  case,  the  inflow  at  the  upper  end  (Wheeling)  and  a 
rating  curve  (a  relation  between  H  and  u  obtained  from  previous 
floods)  at  the  lower  end  (Cincinnati)  were  used.   For  the  second 
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test  case,  the  continuity  conditions  for  the  stages  and  discharges 
at  the  junction  were  used,  in  addition  to  the  boundary  data  at  the 
three  ends.   For  the  third  test  case  the  discharge  rate  at  the 
dams  at  the  two  ends  of  the  reservoir  from  actual  records  were 
used. 

The  two  partial  differential  equations  for  u  and  H  are  a 
hyperbolic  system  possessing  two  real  characteristics.   In  fact, 
the  following  characteristic  form 


Su  , 

St"- 


/    'iA^Sul   igB  r^H   /     /gA^  SH  _  ^^Inl  +  /§  fa   n 


indicates  that  the  two  characteristics  are  given  by  curves  in  the 

dx        fsA  CsA 

X,  t  plane  as  defined  by  -r-p  =  u  ±  J^  =  c_^ .   Thus  J-^  can  be  inter- 
preted as  the  velocity  at  which  a  small  disturbance  in  the  flow 
propagates  downstream  and  upstream  relative  to  the  flow.   The 
hyperbolic  character  of  the  governing  equations  has  an  important 
consequence  in  the  niimerical  scheme  for  solving  the  problem. 
Namely,  for  a  given  space  interval  Ax  for  the  net  points,  the  time 
interval  At  must  be  small  enough  so  that,  j-r   >  c_^_,    with  c_^_,    the 
quantity  defined  above.   In  all  three  test  cases.  Ax  was  chosen 
as  10  miles  and  At  was  taken  as  9  minutes  in  order  to  be  within  the 
necessary  bound. 

Finally,  we  remark  that  the  success  of  the  present  method 
over  the  flood  routing  method  for  some  problems  (for  which  the 
dynamical  equation  of  motion  is  replaced  by  a  kinematic  empirical 
relation  that  is  based  in  part  on  the  idea  that  flood  waves 
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progress  downstream,  thus  backwater  effects  are  ignored)  is  caused 
by  the  fact  that  propagation  of  waves  upstream  as  well  as  down- 
stream has  been  taken  into  account. 
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3-1 
Lecture  3 

Frontal  Motion  in  the  Atmosphere 
J.  J.  Stoker 

The  weather  in  middle  latitudes  is  determined  primarily  by- 
events  associated  with  the  development  and  propagation  of  frontal 
cyclones  between  tropical  and  polar  air  masses  in  the  atmosphere. 
Before  the  advent  of  computers,  mathematical  studies  were  limited 
to  small  perturbations  to  a  quasistationary  front  with  the  aim  of 
finding  solutions  whose  patterns  resemble  those  of  observed 
nascent  cyclones »   As  a  prelude  to  study  the  nonlinear  effects  in 
the  dynamical  equations.  Stoker  (1957)  formulated  a  two-layer 
model  for  the  frontal  motion.   Then,  Kasahara,  Isaacson  and  Stoker 
(1965)  made  the  first  successful  attempt  in  showing  the  evolution 
toward  an  occlusion  from  a  developing  frontal  cyclone. 

Before  describing  the  mathematical  model  and  the  results  of 
numerical  computation,  a  brief  mention  of  some  of  the  meteoro- 
logical observations  about  weather  fronts  should  be  made.   In  the 
middle  latitudes  the  upper  layer  of  warm  air  is  separated  from 
the  lower  layer  of  cold  air  by  a  thin  transition  zone  (a  dis- 
continuity surface)  across  which  both  the  density  and  the  tangen- 
tial velocity  of  the  air  change  abruptly.   The  intersection  of 
this  discontinuity  surface  with  the  ground  forms  the  "front". 
The  front  moves  eastward  as  a  whole.   The  portions  of  the  front 
where  the  cold  air  pushes  the  warm  air  off  the  ground  are  called 
cold  fronts,  and  the  portions  where  the  cold  air  is  receding  from 
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the  warm  air  are  called  warm  fronts.   Since  a  cold  front  usually 
moves  faster  than  its  adjacent  warm  front,  what  is  called  an 
occlusion  may  develop  eventually. 

Figure  1  shows  the  dynamical  system  to  be  studied.   The  cold 
air  at  the  ground  occupies  a  wedge-like  domain  with  warm  air  over 
it.   The  interface  is  given  by  z  =  h(x,y, t),  with  x  in  the  eastward 
direction  and  y  in  the  northward  direction.   The  top  of  the  warm 
layer,  where  atmospheric  pressure  is  essentially  zero,  is  given  by 
z  =  h'(x,y, t).   The  fluid  in  each  layer  is  assumed  to  be  inviscid 
and  incompressible,  with  constant  density  p  and  p'  in  the  cold  and 
warm  layers  respectively,  and  subject  to  gravity.   For  large  scale 
motions  in  meteorology  the  hydrostatic  pressure  law  holds  quite  well, 
and  this  means  that  the  acceleration  of  the  fluid  in  the  vertical 
direction  can  be  neglected.   Hence  the  pressure  in  the  warm  layer  is 
given  by 

p'  (x,y,z,t)  =  p'g(h'  -  z)  , 

and  the  pressure  in  the  cold  air  by  • 

p(x,y,z,t)  =  p'g(h' -  h)  +  pg(h -z)  , 

g  being  the  gravitational  constant.   It  follows  that  the  horizontal 
pressure  gradient  is  independent  of  z.   Since  the  Coriolis  force  due 
to  the  rotation  of  the  earth  may  also  be  regarded  as  independent  of 
z  without  serious  error,  it  follows  that  the  horizontal  acceleration 
of  fluid  particles  is  independent  of  z,  and  hence  that  the  horizon- 
tal components  of  velocities  u,  v  and  u',  v'  will  remain  independent 
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of   z    if   they  were    so   initially,    which   is    assiJined   here.      The    equa- 
tions   of  motion,    in  Eulerian   form,    are 

Su',       ,    Su',       ,    Su',  5h'            „    , 

^  +   ^     ^  +  ^     ^  ^  S  ^  =      ^^      ' 

dt                  dx        dy  ^  dy 
for  the  warm  air,  and 

Su  ,  ,,  ^u  ,  „  Su  ,  „/n'  Sh'  ,  p-p'  Sh 


"5^ 


dx     dy   ^^  p  dx      p   dx^       ' 


Sv  ,    ^v  ,    Sv  ,   /p'  Sh'    p-p'  Shx     „ 
dt     dx     ay   ^^  p  dy      p   dy' 

for  the  cold  air.   Here  f  is  the  Coriolis  parameter,  given  by 
f  =  2(D  sin  (})  with  o)  the  angular  velocity  of  the  earth  and  cj)  a 
reasonable  latitude  for  the  region  under  considerationo   Since  the 
fluids  in  both  layers  are  assumed  to  be  incompressible,  the  equa- 
tions of  continuity  can  be  written  as 

1^  (h'  -  h)  +  J-[u'(h'-  h)]+  ^[v'(h'-  h)]  =0  , 

,J'         The  treatment  can  be  simplified  further  by  assuming  that  u', 
v',  and  h'  do  not  change  in  time.   The  reason  that  the  dynaunics  of 
the  perturbations  in  the  warm  air  relative  to  the  initial  stationary 
state  can  be  neglected  as  a  first  approximation  is  that  the  warm 
layer  can  adjust  itself  through  a  slight  change  in  the  vertical 
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velocity  without  a  substantial  change  in  the  horizontal  velocity. 
Therefore,  when  the  warm  air  is  in  a  stationary  state  given  by 

u'  =  u'  (const.)  ,    v'  =  0  ,    h'  =  -  -  u'y  + const. 

O  o 

the  governing  equations  for  the  cold  air  become 


£. 

-P' 
P 

£ 

-P' 
P 

Sh 

I;  +  u  |x  +  V  |i  +  g  £^  1^  =  f  (^  u'  -  u)  . 

The  interface  between  the  two  layers  is  a  free  surface  in  the  problem. 
Any  point  (x  (t),y  (tj)  on  the  front  moves  with  the  particle  velocity 
on  it.   This  follows  from  the  continuity  equation,  since  at  h  =  0 

the  equation  (-^  "*"  ^  15x  "*"  ^  "Sv^^  "  *^  holds. 

To  solve  the  partial  differential  equations  for  u,  v,  and  h 
a  rectangular  domain  (0<x<X,  0<y<Y)  in  the  xy  plane  is  taken 
(see  Fig.  2).   For  the  boundary  conditions,  v  =  0  is  assumed  for  all 
time  on  the  northern  boundary  (y  =  Y)  and  u,  v,  h  are  assumed  to  be 


periodic  in  x,  viz.   h 


=  h 
x=0 
dx  dy 


etc.   The  boundary  condition  on 
x=X 


the  front  is  that  -ri^-  =   u  and  -j-r^  =  v.   Two  sets  of  initial  condi- 

dt  dt 

tions  were  actually  considered. 
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Case  A.   The  initial  conditions  are: 
u  =  u   (const.)  ,    V  =  0  , 

ye  =  ^2  -  ^1  ^°^  ^  ^c  • 

Case  B.   The  initial  wind  field  is  geostrophic: 
u  =  u^  ,    V  =  0  , 

h  =  ^^  I  -^    (-^  u'  -u  )(Y  -C,  -  C^)  , 
Y-y^,  g  p-p'  ^  p   o    o^^     1    2'    ' 

y   =  C„  -  C-,  cos  -^   X 
"^  c     2    1      X   c 

In  the  computation,  it  is  convenient  to  use  the  moving  coordinates 
^  and  ri  instead  of  x  and  y  as  defined  by 

?  =  X  -u^t  ,    T)  =  y  . 

The  increment  sizes  At  =  10  minutes,  A|  =  76.2  km,  Arj  =  76.2  km  were 
chosen  and  X  =  20A|,  Y  =  20Ar),  C-j_  =  2Ati,  C^  =  9.5Ari,  u  =  lOft/sec, 
-2^  u'  =  50ft/sec,  f  =  10"  sec""^,  g  =  32.15f t/sec^,  g(l  _  £.1.)  =  .6f t/sec^ 
were  used.   The  above  magnitude  of  the  density  discontinuity  corre- 
sponds to  a  temperature  discontinuity  of  5°C.   The  results  of 
numerical  calculation  are  shown  in  Fig.  3  for  Case  A  and  in  Fig.  4 
for  Case  B.   In  Fig,  3c  the  numerators  and  the  denominators  represent 
respectively  the  ^-component  and  the  -Q-component  of  the  velocities 
of  points  on  the  front  relative  to  the  moving  coordinate  system, 
which  has  an  eastward  velocity  of  lOft/sec. 


19 


3-6 


It  is  seen  that  particles  on  the  cold  front  moved  southeast- 
ward and  those  on  the  warm  front  moved  northwestward  on  the  average, 
whereas  both  fronts  themselves  propagated  eastwards.   The  movement 
of  particles  on  the  front  clearly  indicates  the  production  of  a 
cyclonic  circulation  around  a  center  near  the  Junction  of  the  cold 
front  with  the  warm  front.   The  converging  motion  of  the  cold  air 
near  the  circulation  center  will  produce  moist  convections.   This 
implies  that  severe  storms  are  likely  to  be  associated  with  cold 
fronts.   In  conclusion,  the  two  cases  studied  indicate  that  the 
action  of  the  purely  mechanical  forces,  gravity  and  the  Coriolis 
force,  are  primarily  responsible  for  the  gross  features  of  the 
occlusion  process  for  the  weather  fronts. 
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Fig.  1.  {a)  A  stationary  front.  (6)  \'ertical  cros.s  section  of  the  two  layers. 
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Fig.    2.      The    rectangular   domain   of  numerical    integration. 
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Fig.  3«  Case  A.  Height  contours  (at  5000  ft.  intervals 
in  the  cold  air  layer  and  trajectories  of  front 
points  during  the  period  of  8  hours. 
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Fig.  4, 


Case  B.   Height  contours  (at  5000  ft.  Intervals) 
m  the  cold  air  layer  and  trajectories  of  front 
points  during  the  period  of  11  hours. 
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Lecture  4 

Frontal  Motion  in  the  Atmosphere 
E.  Turkel 

The  equations,  derived  in  a  previous  lecture,  for  a  single 
layer  model  of  the  atmosphere  are  to  be  discussed  here.   Those 
equations  are 

U^+UU^  +  VUy  +g(l--^)h^  =  fv 
(1)  v^  +  UV^+VV  +g(l-  £^)h   =  f(u-  ^  u') 


h ,  +  ( hu )  + ( hv )   =0 


where  u  and  v  are  the  velocities  of  the  cold  air  in  the  x  and  y 
directions  respectively,  while  h  is  the  height  of  the  cold  air. 
p  and  p'  are  the  densities  of  the  cold  and  warm  air  masses,  f  is 
the  Coriolis  parameter  and  u'  is  the  speed  in  the  x  direction  of 
the  warm  air. 

It  is  assumed  that  the  cold  air  lies  above  a  region  D  of  the 
X, y  plane  and  that  this  region  is  bounded  on  three  sides  by 
straight  lines  and  on  the  fourth  side  by  a  curve  C(t)  as  illustra- 
ted in  Fig.  1.   The  curve  C(t)  :  (x  (t),y  (t))  is  defined  as  the 
line  along  which  h  =  0.   The  points  of  C  move  in  time  with  the 
velocity  (u„,v_,)  of  the  particles  on  the  front:   Hence  the 
following  conditions  hold  along  C: 
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h  =  0 
(2)  ^   (x^(t))  =  u^(x^,y^,t) 

where  -rr:   is  the  particle  derivative. 

The  quantities  u^  v  and  h  must  be  prescribed  along  the 
other  three  boundaries.   Two  different  situations  are  considered 
for  the  distance  between  the  east  and  west  boundaries.   For 
simplicity  it  is  first  assumed  that  the  dependent  variables  are 
initially  sinusoidal  in  the  x  direction  with  a  wavelength,  1\, 
equal  to  the  distance  between  the  east  and  west  boundaries o   In 
the  second  situation,  the  east-west  boundaries  are  taken  to  be 
three  times  more  distanct  than  in  the  first  situation,  while  the 
Initial  data  are  chosen  to  be  constant  for  0  j<  x  j<  A  and 
2A  <  X  <  3A,  but  in  the  region  A  <  x  <  2A  the  dependent  variables 
are  chosen  to  be  sinusoidal  with  wavelength  A,  so  as  to  lessen  the 
boundary  influence  on  the  occlusion  process.   Along  the  northern 
boiindary  the  north-south  component  of  the  velocity,  v.  Is 
specified;  and  is  not  necessarily  zero  as  was  assumed  in  the 
preceding  lecture. 

A  complete  formulation  of  the  problem  requires  appropriate 
initial  data  for  u,  v  and  h.   It  is  also  necessary  to  specify  the 
initial  position  of  the  curve  C(t)  and  the  values  of  u  and  v  along 
the  front.   In  Case  1,  u  and  v  are  taken  initially  the  same  as  for 
the  constant  steady  state  zonal  solution,  while  h  varies  linearly 
in  y  and  sinusoidally  in  x. 
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The  initial  conditions  are  thus 
Case  1:   (Non-geostrophic ) 


yc  =  ^2-^1  ^°^  ^T^^: 


u  =  u 


V  =  0 


h  =  (y-y2)H  ,       Yq  <y  1^ 


where 


! 


H  =  ^— r-  (£-  u'-u)  ;    £1  u'  >  u  . 

g(l-^)   P  P 

X,  Y  denote  the  lengths  of  the  sides  of  R  in  the  x  and  y  directions 
respectively.  The  contour  lines  of  the  initial  value  of  h,  at  5000 
foot  intervals,  are  plotted  in  Figure  2a. 

As  a  second  case,  a  physically  more  realistic  condition  is 
assiimed,  namely  that  the  wind  field  is  initially  geostrophic. 
Thus,  at  t  =  0  there  is  no  acceleration  in  either  the  x  or  y 
directions,  i.e.  -^   =  -^  =   0.   This  leads  to  the  initial  condi- 
tions. 

Case  2;   (Geostrophic) 

y-yp 

h  =   ^^   (Y-b)H 

f    ay   p 

I      dx 
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where  b  =  C-,  +  Cp  and  H,  Y,  y„  are  as  in  Case  1.   The  contour 
lines  of  the  initial  value  of  h  for  Case  2  is  plotted  in 
Figure  5a. 

Subcases  la  and  2a  denote  the  problems  with  periodic  east- 
west  boundary  conditions  and  v  =  0  at  the  northern  boundary. 
Subcases  lb  and  2b  denote  the  same  problems  but  with  v  specified 
(nonzero)  at  the  north.  Finally,  subcases  Ic  and  2c  concern  the 
problems  with  the  extended  x  domain,  while  v  is  taken  as  zero  at 
the  north. 

Difference  Equations 

Consider  the  rectangle  R  with  sides  L-,,  Lp.   A  dimensionless 
moving  coordinate  system  is  introduced  so  that  in  dimensionless 
units  the  unrefined  mesh  has  length  and  width  equal  to  1  (see 
[3,5])'      By  D.  is  meant  the  connected  set  of  net  points  in  the 
interior  of  D,  the  region  under  the  cold  air.   By  a  regular  point 
is  meant  a  point  in  D  whose  eight  nearest  neighbors  are  all  in  D^. 
All  other  points  in  D   are  called  irregular  points.   At  regular 
points,  two  different  second  order  schemes  are  considered.   The 
first  is  a  one-step  method  that  is  a  nonlinear  generalization  of 
the  Lax-Wendroff  method  [4].   The  second  method  is  a  two-step 
scheme  due  to  Burstein  [2].   At  the  northern  boundary  and  at 
irregular  points  sufficiently  far  from  the  front,  variations  of 
the  Lax-Wendroff  method  were  used.   At  irregular  points  too  close 
to  the  front  the  values  of  u,  v  and  h  were  found  by  interpolation. 

Boundary  Conditions 

The  finite  difference  approximations  at  the  boundaries  are 
considered  next.   The  east  and  west  boundaries  consist  of  regular 
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net  points  upon  using  the  periodicity  condition.   Along  the 
northern  boundary  v  is  given,  while  u  and  h  are  calculated  by 
using  one-sided  differences. 

Along  the  front  the  sound  speed  is  zero  and  so  the 
characteristic  directions  coincide  with  the  particle  path.   Thus, 
u  and  V  cannot  be  found  by  extrapolation  along  the  characteristics 
but  must  be  found  by  solving  the  differential  equations  (2)  for 
u  and  V  along  the  front. 

An  implicit  method,  based  on  the  trapezoidal  rule,  is  used 
to  solve  the  ordinary  differential  equations.   The  accuracy  of 
this  method  is  consistent  with  the  second  order  accuracy  of  the 
Lax-Wendroff  method.   Furthermore,  the  trapezoidal  method  does  not 
require  information  at  previous  time  levels  and  so  presents  no 
difficulties  when  the  points  along  the  front  are  redistributed. 

It  was  found  that  as  the  occlusion  process  continued  the 
individual  particles  on  the  front  converged  toward  the  place  where 
the  warm  and  cold  fronts  join.   Therefore,  the  spacing  between  the 
points  on  the  front  becomes  small  in  comparison  with  the  grid  size. 
It  was  thus  found  necessary,  on  occasion,  to  select  new  equi- 
distant points  along  the  front.   The  values  of  x,  y,  u,  v  at  the 
new  points  on  the  front  are  found  by  interpolation  from  their 
values  at  the  old  points. 

Results 

Figure  2b  shows  subcase  la  after  18  hours.   Evidently  the 
entire  front  progresses  to  the  east  with  the  cold  front  moving 
eastward  faster  than  the  warm  front.   It  is  noticed  that  the  front 
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is  no  longer  a  single-valued  function  of  the  x  coordinate  and 
that  the  front  is  curling  counterclockwise.   The  development  of 
this  asymmetry  suggests  the  occlusion  process  of  frontal  cyclones. 

Subcase  la  assumes  that  no  air  enters  from  the  north.   The 
problem  was  then  solved  for  subcase  lb  where  cold  air  enters  the 
region  from  the  north  simulating  a  polar  wind.   As  in  Alterman 
and  Isaacson  [1],  the  velocity  v  was  assumed  given  by 

>p[cos  (I^)-  1]    ,         t  <  2T  , 


(3)  V(t,Y,x) 


0  ,    t  >  2T  . 


Figure  3  shows  the  result  after  18  hours  with  T  =  7  hours  and 
V  =  10  ft/sec.   As  is  natural,  it  is  seen  that  the  polar  air 
forces  the  front  southwards   With  cold  air  entering  from  the 
north,  the  front  becomes  steeper  and  the  cyclonic  activity 
increases. 

The  periodic  east-west  boundary  conditions  are  not  realistic 
for  a  single  wave  of  period  X.   Thus,  in  subcase  Ic  a  larger 
region  was  introduced  to  simulate  an  infinite  domain  in  the  x 
direction.   That  is,  X  was  chosen  as  6o  dimensionless  units; 
while  the  initial  data  was  chosen  to  be:   a  single  wave  of  period 
^  over  the  center  twenty  units  and  extended  to  the  right  and  left 
as  having  values  independent  of  x,  see  Figure  4a.   As  seen  in 
Figure  4b  the  main  occlusion  process  is  almost  unchanged  by  the 
new  boiindary  conditions  after  18  hours.   In  Figure  4c  we  plot  the 
contour  lines  of  h,  for  the  midportion   (20  jf  x  ^  4o),  to  facili- 
tate their  comparison  with  the  corresponding  lines  in  Figure  2b. 
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Note  that  in  the  regions  away  from  the  center  of  the  occlusion 
process,  the  front  remains  closer  to  its  initial  state  and  even 
moves  slightly  southward.   The  propagation  of  disturbances  is 
faster  towards  the  north  since  the  sound  speed  is  proportional  to 
the  square  root  of  the  height  of  the  cold  air.   This  height 
increases  approximately  linearly  towards  the  north. 

In  Case  2  the  initial  velocities  and  height  field  distribu- 
tion were  assumed  to  be  geostrophic.   The  initial  shape  of  the 
front  is  again  sinusoidal  with  the  height  contour  lines  flattening 
out  towards  the  north  as  is  seen  in  Figure  5a.   Figure  5b  shows 
the  position  of  the  front  together  with  the  height  contour  lines 
after  24  hours.   As  before  the  entire  frontal  system  progresses 
eastward.   It  takes  more  time  for  the  occlusion  process  to  occur 
since  the  initial  accelerations  are  zero.   However,  after  a 
complete  day  the  occlusion  process  is  apparent  and  large  cyclonic 
motions  have  developed. 

Subcase  2b  considers  the  effect  of  a  polar  wind  entering 
from  the  north  as  given  by  equation  (3).   Figure  6  shows  the 
contour  levels  for  the  height  field  for  this  case  after  24 
hours.   As  before  the  boundary  conditions  were  then  changed  with 
the  east-west  boundaries  taken  as  three  times  further  apart  and 
with  no  wind  entering  from  the  north.   In  Figure  7a  is  shown  the 
initial  configuration  for  subcase  2c.   In  Figure  7b  it  is  seen 
that  the  occlusion  process  is  not  greatly  affected,  after  24  hours, 
by  this  change  in  boundary  conditions.   Since  the  height  levels 
are  relatively  independent  of  x  toward  the  north  the  wave  motion 
is  not  as  well  developed  as  it  was  in  Case  1.   Figure  7c  shows 
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the  center  portion  of  Figure  7b  to  facilitate  comparisons  with 
Figure  '^h. 

Conclusions 

The  general  motion  of  cold  fronts  in  the  atmosphere  can  be 
qualitatively  described  by  a  system  of  nonlinear  equations  for 
the  motion  of  a  single  layer  of  fluid.   These  equations  include 
gravitational  and  Coriolis  forces  but  ignore  thermodynamic 
variables.   By  changing  the  initial  and  boundary  conditions,  the 
shape  of  the  front,  the  cyclonic  pattern  of  the  air  and  the  time 
variation  can  all  be  varied.   Within  one  day  the  cyclonic  pattern 
is  clearly  developed  for  realistic  initial  conditions.   Thus,  this 
model  indicates  that  the  action  of  purely  mechanical  forces,  i.e. 
gravity  and  the  Coriolis  effect,  are  primarily  responsible  for  the 
gross  features  of  the  occlusion  process. 
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Fig.  2a:   Initial  height  contour  pattern,  in  5  thousand  foot  intervals, 
for  case  la  with  no  air  entering  from  the  north. 
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Figo  4c:   Central  portion  of  Fig..  4b  to  facilitate  comparison 
with  2"b . 
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Fig.  5b:   Height  contour  pattern  for  case  2a  after  24  hours. 
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Lecture  5 

Mountain  Winds 
E.  Isaacson 

The  occurrence  of  dangerous,  high  speed  winds  on  the  lee 
slope  of  a  mountain  (e.g.  the  east  slope  of  the  Rocky  Mountains) 
is  still  unpredictable.   To  understand  such  strong  mountain  winds, 
Houghton  and  Kasahara  ( I968 )  studied  a  one-layer  atmospheric 
model  over  a  symmetrical  mountain  to  detemine  the  steady  asymptotic 
states  of  the  mountain  flows  arising  from  a  family  of  problems  in 
which  the  west  to  east  wind  is  perpendicular  to  the  mountain  ridge 
and  instantaneously  constant  at  the  initial  time.   They  found  that 
if  this  initial  velocity  is  large  enough  (but  not  extremely  large), 
then  the  flow  is  asymmetrical  and  the  asymptotic  state  contains  a 
hydraulic  jump.   Such  a  jump  either  remains  stationary  on  the  lee 
side  of  the  mountain  or  moves  downwind  with  a  small  constant 
velocity  depending  on  the  magnitude  of  the  initial  velocity.   This 
one-layer  model  gives  an  excellent  explanation  of  the  hydraulic 
jump,  often  seen  as  white  turbulence,  on  the  downstream  side  of  a 
large  boulder  immersed  in  a  shallow,  rapid  stream.   In  order  to 
simulate  nature  more  closely,  Houghton  and  Isaacson  (1969)  studied 
a  two-layer  model  of  the  atmosphere.   Their  aim  was  to  find  the 
asymptotic  states  of  the  asymmetrical  flows  that  may  arise  and  to 
compare  these  with  observed  flows. 

Based  on  observations  by  Harrison,  the  atmosphere  between 
the  surface  of  the  earth  and  the  tropopause  is  assumed  to  consist 
of  two  layers,  separated  by  a  very  thin  mid-tropospheric  interface. 
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Each  layer  is  treated  as  an  incompressible  fluid  subject  to  the 
hydrostatic  pressure  law.   The  lower  layer,  between  the  ground, 
z  =  H(x),  and  the  mid-troposphere  interface,  has  a  depth  i{x)    and 
an  average  horizontal  flow  velocity  u(x)  with  a  constant  density  p. 
Whereas  the  upper  layer,  between  the  mid-tropospheric  interface 
and  the  tropopause,  has  a  depth  (l)'(x)  and  an  average  horizontal 
flow  velocity  u' (x)  with  a  constant  density  rp  (r  <  l).   The 
stratosphere  on  top  of  the  tropopause  is  assumed  to  be  at  rest  with 
a  horizontal  upper  surface.   This  passive  stratospheric  layer  has 
a  constant  density  sp  (s  <  l).   The  governing  equations  for  (f),  u, 
(f)',  and  u'  are 

Ie  *  +  |j  (*u)  =  0  , 

1;  u.u  |j  u.g(r-B)  1^  V+  g(l-s)  1^  *  =  -g(l-s)  i  , 

^  u.  .u.  |j  u.  .g(i- 1)  |j  *..  g(i-|)  |L  i  ,  .g(i.  1)  as  , 

g  being  the  acceleration  of  gravity. 

Here,  for  the  case  of  a  symmetrical  mountain,  the  nature  of 
the  asymptotic  states  (as  t  -»  co  ),  that  arise  from  initially 
constant  velocities  and  horizontal  atmospheric  interfaces,  is  more 
complex  than  for  the  case  of  the  single  layer  flow.   In  fact, 
except  for  the  case  of  quite  low  initial  velocity,  in  which  the 
flow  is  again  symmetrical,  there  is  a  regime  of  steady  asymmetric 
flow  over  the  mountain,  only  for  a  narrow  neighboring  range  of 
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higher  initial  velocities  (region  B  In  [2]).   If  the  initial 
velocity  Is  chosen  to  be  still  somewhat  higher  (region  B'  in  [2]), 
then  the  flow  over  the  mountain  falls  to  become  steady,  but  instead 
behaves  roughly  as  if  periodic  waves  are  created  over  the  mountain 
and  sent  downstream.    Nevertheless,  both  this  steady  asymmetric 
flow,  with  its  hydraulic  jumps  on  the  lee  side,  and  the  unsteady 
flow,  occur  for  physically  realistic  atmospheric  velocities. 

With  this  study  of  flow  over  a  symmetrical  mountain  completed, 
a  calculation  of  the  flow  over  a  realistic  profile  of  the  Rocky 
Mountains  was  undertaken.   The  ground  is  assumed  to  have  the 
average  profile  of  a  ten  mile  wide  strip,  running  east-west  through 
Colorado  at  a  latitude  between  Denver  and  Boulder,  Colorado.   If 
initially  each  of  the  two  layers  has  a  horizontal  upper  surface 
and  a  uniform  eastward  flow  velocity  as  indicated  in  Fig.  1,  the 
state  after  ^.IJ   hours  is  shown  in  Fig.  2  for  the  heights  of  the 
upper  surfaces  of  the  two  layers  and  the  flow  velocities  in  the  two 
layers.   The  solution  shows  a  severe  decrease  in  depth  of  the  lower 
layer  on  the  lee  side  of  the  mountain,  and  a  correspondingly  high 
wind  velocity,  u,  in  the  Boulder  area.   In  fact,  u  is  100^  larger 
than  u';  and  u  is  close  to  the  125  niph  winds  observed  at  Boulder 
in  [5]  .   In  [2],  a  comparison  is  made  of  another  calculated  flow, 
with  carefully  made  observations  of  the  atmosphere  on  Feb.  20,  I968. 


Subsequent  work  by  Isaacson  and  Zwas  [4]  shows  that  this  unsteady, 
almost  periodic  behavior  was  caused  by  the  numerical  scheme.   In 
fact,  the  flow  for  such  data  (called  B'  in  [2])  does  become  steady, 
but  it  contains  a  secondary  weak  shock  on  the  lee  side  of  the  crest. 
The  numerical  schemes  to  carry  out  this  calculation  are  studied 
in  [4]. 
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It  is  pointed  out  that  this  model  is  easy  to  use  for  local  wind 
forecasting  purposes.   Further  high  wind  observations  are  needed 
to  justify  more  fully  the  use  of  this  two-layer  model,   which 
predicts  the  occurrence  of  a  high  shear  in  the  wind. 
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Lecture  6 

Atmospheric  Predictability 
(an  introduction  to  large-scale  numerical  meteorology) 

Richard  C.  J,  Some rvi lie 

Introduction 

Suppose  that  one  possessed  perfect  knowledge  of  the  physical 
laws  which  describe  the  behavior  of  the  earth's  atmosphere  and  of 
the  appropriate  mathematical  expressions  of  these  laws.  Suppose 
too  that  the  numerical  methods  smd  computational  machinery 
necessary  for  an  essentially  errorless  integration  of  the  re- 
sulting initial-  and  boundary- value  problem  were  available.  Also 
suppose  that  one  could  obtain  all  the  observational  data  needed 
to  define  the  initial  and  boundary  conditions.  With  such  perfect 
models,  methods,  and  data,  for  how  long  a  time  could  one 
"accurately"  (defined  in  some  appropriate  sense)  forecast  the 
weather? 

Alternatively,  for  the  imperfect  models,  methods,  and  data 
actually  available  today,  in  what  respects  are  forecasts  deficient, 
and  for  how  long  a  time  can  forecasts  presently  be  said  to  be 
accurate? 

Finally,  to  produce  an  accurate  forecast  for  some  arbitrary 
length  of  time,  intermediate  between  the  presently  practicable 
and  the  ultimately  possible  lengths  of  time,  what  sorts  of  models, 
methods,  and  data  would  be  required? 
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Atmospheric  predictability  is  the  subject  which  is  concerned 
with  such  questions.   It  is  a  subject  of  practical  importance  as 
well  as  theoretical  interest,  because  while  the  economic  benefits 
of  improved  weather  forecasts  are  immense,  so  are  the  costs  of 
super-computers,  satellites,  etc.   In  the  following  sections,  we 
shall  first  describe  the  problem  a  little  more  concretely  and  then 
briefly  survey  several  separate  lines  of  attack  on  it.   The  survey 
is  by  no  means  exhaustive,  but  the  bibliographies  contained  in  the 
references  provide  an  extensive  set  of  entry  points  to  the 
literature, 

Basic  Equations 

The  appropriate  fundamental  laws  are  thought  to  be 
representable  by  the  follov;ing  system: 

^V  =  -  SnxV-  aVp  +g  +F 

■g^  a  =  ay.V 


^T  =  (1-7)TV.V+| 


V 


pa  =  RT  . 


These  are,  respectively,  Newton's  second  law,  the  statement  that 


mass  is  conserved,  the  first  law  of  thermodynamics,  and  an  equation 
of  state.  The  individual  time  derivative  -^  means  (-j^+V*^)*  Z  ^^ 
the  velocity  relative  to  the  earth  rotating  at  angular  velocity  O, 
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a  is  the  reciprocal  of  density,  p  is  pressure,  g  is  the  apparent 

gravitational  acceleration,  R  is  the  gas  constant,  y   is  the  ratio 

c 
of  specific  heats   {->£),  and  T  is  temperature.   The  frictional 

*^v 
force  P  and  net  heating  Q  per  unit  mass  are  very  complicated  and 

imperfectly  known  functions  of  other  variables  (cf.  Lorenz,  I967). 

Since  World  War  II,  much  effort  has  been  devoted  to  the 
numerical  solution  of  these  equations,  with  a  great  variety  of 
domains,  boundary  conditions,  numerical  methods,  initial  data, 
assumptions  about  F  and  Q,  etc,  (e.g.,  Thompson,  1959;  Gavrilin, 
1965;  Phillips,  1970).   Today,  such  models,  run  on  a  routine 
operational  basis,  comprise  the  core  of  the  process  by  which 
regular  weather  forecasts  are  produced  in  the  technologically 
advanced  countries.   The  more  complicated  models,  which  cannot 
meet  the  operational  real-time  requirement,  are  used  for  research 
purposes.   Quite  aside  from  the  fact  that  the  archive  of 
meteorological  observations  constitutes  what  must  be  among  the 
largest  data  sets  in  all  science,  the  sheer  computational  magni- 
tude of  these  efforts  is  awesome  by  conventional  standards.   The 
largest  models  (general  circulation  models)  saturate  the  memories 
of  the  largest  computers  and  may  require  several  hundred  hours  of 
machine  time  for  a  single  experiment.   The  verisimilitude  of  the 
results  is  sometimes  purchased  at  a  steep  price,  not  only  in 
dollars  but  also  in  loss  of  insight  and  intelligibility  —  the 
model,  like  the  atmosphere  itself,  may  be  so  complicated  that 
specific  effects  cannot  easily  be  attributed  to  specific  causes. 
Before  examining  the  results  of  the  large  models,  therefore,  we 
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shall  first  survey  three  somewhat  simpler  approaches  to  the 
predictability  problem. 

Analogues 

We  shall  first  briefly  mention  some  fluid  systems  designed 
to  resemble  the  atmosphere  in  certain  important  respects  but  to 
differ  from  it  in  size,  simplicity,  and  susceptibility  to  observa- 
tion and  control.   The  most  extensively  studied  of  such  systems 
are  the  "rotating  annulus"  experiments  (Fultz,  et  al.,  1959), 
which  consist  of  two  concentric  right  circular  cylinders  mounted 
on  a  turntable  and  rotated  together  at  a  constant  rate  about  their 
common  (vertical)  axis.   If  the  cylinders  are  maintained  at 
different  constant  temperatures,  water  contained  in  the  space 
between  them  has  been  found  to  exhibit  many  of  the  features  of 
large-scale  atmospheric  flow,  for  certain  ranges  of  values  of  the 
external  parameters.   Many  of  these  features  (which  include  fronts, 
jet  streams,  and  "planetary"  waves)  have  been  shovm  convincingly 
to  bear  a  more  than  superficial  resemblance  to  their  atmospheric 
counterparts.   These  features,  in  both  annulus  and  atmosphere, 
are  due  to  the  horizontal  temperature  contrast  and  the  dominant 
constraint  of  rotation  (Greenspan,  I968).   Several  qualitatively 
different  regimes  of  flow  can  be  produced  in  the  annulus  by 
varying  the  external  parsimeters.   These  range  from  regimes  of 
infinite  predictability  (steady  flows  and  flows  which  appear  to 
vary  perfectly  periodically  in  time)  to  irregular  and  turbulent 
regimes  which  are  often  supposed  to  correspond  most  closely  to  the 


55 


6-5 


atmosphere.   The  annulus  experiments  have  provided  insight  into 
the  dynamics  of  the  atmosphere,  but  the  analogy  is  too  crude  to 
supply  detailed  quantitative  answers  to  questions  of  predicta- 
bility. 

A  very  different  kind  of  analogue  may  be  sought  in  the 
meteorological  records  themselves.   The  idea  is  to  find  two 
actual  observed  states  of  the  atmosphere,  widely  separated  in 
time,  which  resemble  each  other  so  closely  that  the  differences 
between  them  are  typical  of  the  small  errors  which  might  be  due 
to  an  imperfect  observing  system.   Following  the  subsequent 
history  of  these  two  states  might  then  reveal  the  evolution  of 
the  "error,"  an  estimate  of  the  doubling  time  of  such  errors,  and 
ultimately  the  predictability  time  (after  which  the  two  histories 
are  essentially  uncorrelated,  the  states  of  one  differing  from 
those  at  the  corresponding  times  of  the  other  by  the  same  amount 
as  would  two  randomly  chosen  states). 

A  search  of  the  existing  aerological  records  by  Lorenz 
failed  to  produce  any  such  pairs  of  states  but  did  produce  pairs 
of  states  which  differed  from  one  another  by  amounts  which  might 
be  typical  of  errors  larger  than  errors  of  observation.   The 
observed  growth  rates  of  these  "errors"  could  be  extrapolated  to 
yield  estimates  of  the  growth  rate  of  smaller  errors.   The  result 
(a  doubling  time  of  under  three  days)  is  not  incompatible  with 
results  obtained  by  other  methods,  which  we  shall  now  describe. 
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Simple  theoretical  models 

The  mechanism  implicit  in  the  previous  discussion  (a  limita- 
tion on  predictability  due  to  the  instability,  growth,  and 
eventual  dominance  of  initially  small  observational  errors  in  a 
turbulent  atmosphere)  has  been  further  explored  in  several 
theoretical  studies.  We  shall  briefly  describe  that  of  Lorenz 
(1969).   He  considers  a  flow  governed  by  the  inviscid  two- 
dimensional  barotropic  vorticity  equation,  which  in  terms  of  a 
stream  function  ^  is 

This  equation,  while  much  simpler  than  the  system  described 
earlier,  served  as  the  successful  basis  for  early  short-term 
nximerical  weather  forecasts.   From  it  Lorenz  derives  expressions 
for  ensemble  averages  of  quantities  identifiable  as  the  "error 
energy"  of  separate  scales  of  motion.   After  making  several 
statistical  assumptions  and  simplifications,  a  system  of  ordinary 
differential  equations  (in  time)  results,  the  coefficients  of 
which  may  be  evaluated  from  estimates  of  the  atmospheric  energy 
spectrum.   Numerical  integration  of  this  system  is  easy,  and 
estimates  of  predictability  times  are  obtained  by  noting  the  times 
at  which  the  error  energy  of  a  particular  scale  of  motion  reaches 
the  amplitude  of  the  total  energy  of  that  scale.   With  the  chosen 
values  of  the  coefficients,  motions  on  the  scale  of  cumulus  clouds 
can  be  predicted  for  about  an  hour,  "synoptic-scale"  motions 
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(those  of  typical  weather  map  features,  several  hundred  miles  in 
size)  for  a  few  days,  and  the  largest  (planetary)  scales  for 
"a  few  weeks."  The  predictability  times,  however,  depend 
crucially  on  the  estimates  of  the  atmospheric  energy  spectrum, 

and  may  be  much  longer  if  energy  per  unit  wavenumber  k  falls  off 

-S/5 
faster  than  the  k  '" -^   Kolmogoroff  spectrum  asstimed  for  short 

wavelengths.   Neither  observations  nor  present  theories  of  quasi- 
two-dimensional  turbulence  are  sufficient  to  determine  the 
spectrum  shape  with  any  certainty,  but  it  is  almost  surely 
steeper  than  the  Kolmogoroff  spectrum  (Charney,  I969).   To 
escape  the  necessity  of  such  assiimptions  and  simplifications, 
one  must  resort  to  the  numerical  models.   In  this  subject,  simple 
approaches  provide  insight,  but  realistic  quantitative  answers 
seem  to  come  only  from  less  tractable  procedures. 

Numerical  models 

As  an  example,  we  cite  the  results  of  Smagorinsky  (1969), 
who  describes  his  model  as  follows: 

"Briefly,  the  model  is  governed  by  the  primitive  equations 
and  has  9  vertical  levels.   The  domain  is  hemispheric  and  the 
computational  grid  is  mapped  stereographically.   The  model  is 
'moist'  with  a  condensation  criterion  of  80^  relative  humidity; 
the  parameterized  convection  scheme  is  a  moist  adiabatic  tempera- 
ture adjustment.   Long  and  short  wave  radiation  are  accoxinted  for, 
with  water  vapor,  carbon  dioxide  and  ozone  as  absorbers  of  radiant 
energy  which  at  most  are  functions  of  height  and  latitude  and 
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constant  with  time.   Large  scale  mountains  and  the  thermal  con- 
sequences of  land-sea  contrast  are  accounted  for.  The  surface    .^ 
drag  coefficient  is  constant  irrespective  of  land  and  sea;  the 
availability  of  water  for  evaporation  is  0.5  over  land  and  1.0    ,.. 
over  sea;  the  temperature  specification  is  different  over  land-ice 
and  sea  ice;  the  snow  line  is  constant  with  time.  The  effective 
Karman  constajit  for  the  internal  non- linear  lateral  viscosity 
is  0.4." 

Using  this  model,  with  horizontal  resolutions  as  fine  as  4o 
grid  points  between  pole  and  equator  (a  grid  size  of  320  km  at  the 
pole  and  l60  km  at  the  equator),  pairs  of  3-week  integrations  were 
carried  out.   One  run  of  a  pair  used  initial  data  from  observa- 
tions for  a  particular  day  and  time,  and  the  other  used  this  data 
set  perturbed  by  a  rajidom  temperature  disturbance  of  0.5  deg  C 
amplitude  (standard  deviation).   Thus  an  attempt  was  made  to 
simulate  the  situation  which  Lorenz  sought  unsuccessfully  in  the 
actual  atmosphere.   The  "error"  (instantaneous  standard  deviation 
temperature  difference  between  the  states  of  the  pair  of  runs,  at 
corresponding  times)  first  decreased  from  0.5  to  0.2  due  to  an 
internal  adjustment  with  the  undisturbed  wind  field.   Thereafter, 
the  error  increased  exponentially  for  about  a  week,  with  a 
doubling  time  of  about  2i  days,  and  then  grew  more  slowly.  After 
21  days  it  was  still  appreciably  smaller  than  the  difference 
between  two  randomly  selected  states  at  the  same  season. 

Those  who  perform  such  niamerical  experiments  interpret  them 
to  imply  that  the  ultimate  limit  of  deterministic  predictability 
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may  be  at  least  5  weeks.   The  present-day  practical  limit,  deduced 
by  similarly  comparing  such  forecasts  with  the  actual  weather,  is 
about  half  that  time.   It  is  hoped  that  improvements  in  models 
and  data  will  help  to  close  this  gap,  and  that  increases  in 
computer  capability  will  eventually  allow  routine  operational 
predictions  for  such  extended  ranges. 
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Lectures  7  and  8 

Geostrophic  Vortex  Motion  and  Applications 
G.  K.  Morikawa 

The  large-scale,  nearly  horizontal  motions  of  the  earth's 
atmosphere  are  approximately  in  hydrostatic  and  geostrophic 
balance.   That  is,  1)  in  the  vertical  direction,  the  fluid 
acceleration  is  small  compared  to  both  the  acceleration  of  gravity 
and  the  vertical  pressure  gradient,  viz.,  there  is  approximate 
hydrostatic  equilibrium,  and  ?)  in  the  horizontal  direction,  the 
fluid  acceleration  is  small  compared  to  both  the  Coriolis  accelera- 
tion due  to  the  earth's  rotation  and  the  horizontal  pressure  gradi- 
ent, viz.,  the  motion  is  approximately  geostrophic.   The  horizontal 
balance  is  valid  primarily  in  middle  latitudes,  away  from  the 
equator  where,  strictly,  there  is  no  Coriolis  effect.   If  we 
incorporate  these  meteorological  approximations  in  a  proper  way 
into  the  equations  of  motion  of  an  inviscid,  adiabatic  atmosphere, 
we  can  derive  simpler,  more  manageable  equations  which  succinctly 
describe  these  large-scale  motions;  this  turns  out  to  be  an 
asymptotic  description  for  large  time  of  the  original  equations. 
This  approximate  theory  of  geostrophic  vortex  motion  can  be 
applied  to  problems  related  to  the  study  of  short  range  weather 
forecasting^   ,  valid  for  predictions  over  a  period  of  a  few  days. 
Tvro  such  problems  will  be  considered  here.  We  remark  that  this 
theory  of  kinematic  motion  (dynamically  balanced)  may  be  appli- 
cable to  other  physical  Phenomena  over  a  large  range  of  length 
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scales,  i.e.,  from  atomic  scale  motions  such  as  those  that  occur 

(2  ) 

in   superconductivity  ^  '     to  astrophysical  motions  such  as  inter- 
acting galactic  nebulae. 

Derivation  of  the  Geostrophic  Vorticity  Equation ^-^ 

In  order  to  simplify  the  derivation  of  the  geostrophic 
conservation  equation,  we  consider  a  simpler  model,  that  of  a 
homogeneous,  incompressible  atmosphere  as  our  original  model;  an 

analogous  conservation  equation  results  from  the  geostrophic 

(h  ) 
approximation  of  the  adiabatic  atmosphere''   .   In  the  homogeneous 

atmosphere  model,  the  hydrostatic  approximation  has  already  been 

incorporated;  these  shallow  water  (or  long  wave)  equations  are: 

(1)     —  +   h(u^+Vy)  =  0  (continuity) 

(2.1)  •57-+  ghx-fv  =  0  (x-momentum) 

(2.2)  ■gr  ■•"  gh„  +  f u  =  0  (y-momentum) 

where  h  is  the  depth  of  the  homogeneous  atmosphere,  (u, v)  is  the 

horizontal  velocity,  g  is  the  acceleration  of  gravity, 

f  =20  sin  <j)  (Og  =  earth's  rotational  velocity  and  ^   =  latitude 

angle),  and  d(  )/dt  =  (  )^.  +  u(  )^+v(  )   (subscripts  are  partial 

L      X      y 

derivatives).   For  describing  meteorological  motions,  we  replace 
the  continuity  equation  (l)  by  the  potential  vorticity  equation 
obtained  by  combining  (l)  with  the  curl  of  (2): 
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(3)     ^  ^^TT^   ^   °  (potential  vorticity) 

where  ^  =  (v  -  u  )  is  the  vertical  component  of  vorticity.   Equa- 
tions (2)  and  (3)  [or  (l)  and  (2)]  are  a  system  of  hyperbolic 
equations  for  (h, u, v).   In  addition,  the  pressure,  p,  and  vertical 
velocity,  w,  are  given  by 

(4.1)  p  =  gp(h-z)  (hydrostatic  equilibrium) 
and 

(4.2)  w  =  -z(u^-(-Vy)  . 

We  incorporate  the  geostrophic  approximation  into  (2)  and 
(3)  [or  (l)  and  (2)]  by  the  following  formal  procedure:   l)  the 
time,  t,  is  scaled  by  a  small  parameter,  e,  where  e  <<  1;  and 
2)  the  solution  (h, u, v)  is  expressed  as  a  power  series  expansion, 
perturbed  on  the  atmosphere  at  rest: 

(5)     (h,u,v)  =  {h,0,O)  +  e{h^^\n^^lv^^h  +  oU^) 


where  h  =  constant,  the  height  of  the  homogeneous  atmosphere  at 
rest.   The  first  order  approximation  for  (h^   ,u^   ,v^  ')  yields 


,(1) 

It 

,(1)  _  r.M)    = 


(6.1)   lg_[(c(^)  +  fh^^))/h^]  =  0  ,    T  =  St  ,    f  =  const. 


(6.2)  gh^l^  -  fv^-^^  =  0  [ 

(6.3)  gh^^^  +  fu^^^  =  0  j 


(geostrophic  balance) 
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where  d^^^  )/dT  ^  (  )  +u^^^  )^ +v^^^  ).     We  call  (6.1)  the 
"geostrophic  conservation  equation".   The  system  of  equations  (6) 
can  be  written  in  a  neater  form  by  replacing  h^  ^  by  ^  =  gh^  '/f: 

d(l)      2 

(7.1)  ^  {/^-K^h!^   =  0 

(7.2)  (u^^^v^^^)  =  (-^y,^^) 

where  a(  )  =  (  )   +(  )   *  is  the  two  dimensional  Laplacian,  and 

XX     yy 

K^  s  f^/gh  .   The  geostrophic  equation  (7.2)  shows  that  f   is  the 
stream  function,  viz.,  the  flow  is  divergence-free  to  this  order 
of  approximation;  and  the  first  order  vertical  velocity  w^  '  =  0 
from  (4.2).   If  we  compare  the  asymptotic  equation  (6.1)  [or 
(7.1)]  with  (3),  we  see  that  the  essential  consequence  of  the 
approximation  procedure  is  that  the  potential  vorticity,  (^+f)/h, 
has  been  linearized  to  yield  (^^ ''"^  +  fh^"""  Vh^)  -  U  -  ic    hl^i    however, 
(7.1)  retains  the  nonlinearity  of  the  first  order  material  deriva- 
tive operator,  d^  '(  )/dT.   For  the  linearized  potential  vorticity 
operator,  we  can  obtain  singular  solutions  which  we  call  recti- 
linear  geostrophic  vortices;  such  a  single  vortex  at  r  =  r   is  the 
solution  of 

(8.1)  (A-/c^)^Q  =  5(|r-r^|) 

where  |r-r^I  =  [  (x-x^)^  +  (y-y^)^]"""^^.   The  solution  of  (8.1)  is 


->  -> 


(8.2)  i,^   =  .^Kj.Ir-rJ) 
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where  'y/2iT   is  the  vortex  strength  in  v/hich  v;e  have  chosen  the  sign 
of  7  so  that  positive  7  corresponds  to  a  counter-clockwise 
rotating  vortex  (cyclonic,  in  meteorological  terminology),  and  K 

is  a  Bessel  function.   The  vortex  solution  (8.2)  alone  is  a 
stationary  solution  of  (7);  in  fact,  any  axially  symmetric  solu- 

tion  of  (7)  is  stationary.   Since  {ts  -  ic    )ip     =0  everywhere  except 

->• 
at  r  ,  the  vorticity  distribution  of  a  geostrophic  vortex  is 

(8.3)  ^^^^  =  A'^^  =  /c\  =  -  If-  ^o^'^'^-^oI) 

and  t,^    '   and  ip     have  the  same  sign.   With  respect  to  the  vortex 

-> 

center,  r  ,  the  tangential  velocity  distribution  is 

(a.4)       ^..gK;w;.;j)=gK,ui;-?j)  . 

Near  the  origin, (r  — >  r  ), 


^  or     I"*  ^ 


-  *^  — 

'  r-r. 


o ' 


which  behaves  like  the  classical  logarithmic  vortex.   At  large 

distances  from  r   (  r-r   — >  00  )  , 
o      o 

(8.6)        ^^  (dr-r'j)^/^  exp  [-(^r-rj)] 


which  dies  out  faster  than  the  classical  vortex.   Equations  (8) 
show  that  the  geostrophic  vortex  has  the  qualitative  features  of 
a  closed  low  (or  high)  pressure  system,  such  as  a  hurricane.   We 
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mention  here  that  with  an  adiabatic  atmosphere  we  can  obtain  a 

(4) 
geostrophic  vortex  which  has  vertical  structure. 


Representation  of  a  Hurricane  by  a  Rectilinear  Geostrophic 
Vortex*"^ -^  , 

We  can  use  (7)  in  the  problem  of  short-range  weather  predic- 
tion in  middle  latitudes  of  large-scale  atmospheric  motions,  far 
enough  above  the  ground  where  real  gas  effects,  such  as  viscosity, 
can  be  ignored.   Ordinarily,  this  involves  numerical  calculations 
on  electronic  computers  since  (7.1)  is  a  third-order  nonlinear 
partial  differential  equation.   The  availability  of  a  vortex  solu- 
tion (8)  facilitates  such  calculations  in  the  problem  of  following 
hurricane  trajectories.  We  write  ?//  as  the  sum  of  two  terms 


(9) 


V/(x,y,T)  =  Tp^{i^\'r-r^\) +  i/-^{x,Y,x) 


where  ^   is  given  by  (8.2)  with  r  =r  (t)  and  represents  the  hurri- 
cane  centered  at  r   and  ijj-.    is  the  remaining  part  of  the  flow  field. 
If  we  put  (9)  in  (7),  we  obtain  (as  in  the  classical  vortex  motion 
case ) 


(10.1) 

and 

(I0o2) 


d(l)      2 


"ar 


f 


dx   dy 
o  ^_o 

dx  '  dx 


1 


(xQ,y^) 


where  d 


(1), 


)/dT=^. 


u 


(1)  S( 


"ST 


•+  V 


(1)  hi  )     (1) 


'^T' 


u 


If(^o^^i)- 
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f  1 )    S 
and  v^  '^  =  -^  {tIj    +Tp-.).      The  point  [x  (T),y  (t)]  is  taken  as  the 

position  of  the  hurricane  eye  at  time,  t.   Equations  (lO)  are  a 
coupled  system  of  a  partial  differential  equation  for  f,    and  two 
ordinary  differential  equations  for  (xQ,y  ).   Given  the  initial 
conditions  V-i(^*y>0)>  ^o^^^  ^"^  ^o^^^  (and  boundary  conditions) 
we  obtain  ^.(x,  y,  t)  and  [x^(t),  ^^^(t)]  at  any  later  time   t, 
usually  by  numerical  computation.   We  show  the  results  of  calcula- 
tions  for  Hurricane  Betsy^   ,  13-19  August,  1956.   During  this 
time  the  hurricane  was  in  a  relatively  mature  state  and  the 
trajectory  stayed  over  the  Atlantic  Ocean  without  landfall.   The 
initial  data  (originally  from  the  U.S.  Weather  Bureau)  was  taken 
at  the  500  mb .  pressure  level  (ground  level  is  approximately 
1000  mb.)  and  prepared  in  the  Meteorology  Department,  University 
of  Chicago.   The  height  data  was  available  on  a  22x  22  mesh,  with 
300  km.  mesh  width  at  12  hour  intervals.   The  original  runs  were 
made  with  one  hour  time  steps  for  48  hours  on  the  UNIVAC  computer 

and  later  rerun  at  one-half  hour  time  steps  on  the  IBM  7C4  and 

(7)  (d) 

7090.   The  finite  difference  scheme^  '  and  computer  program''  ' 

were  developed  by  Eugene  Isaacson  and  David  Levine.   Using  the 

available  initial  data  at  each  12  hour  interval,  two  runs  were 

made:   l)  one  run  with  7=0  (non- interacting  vortex)  and  2)  one 

run  with  y   and  k   chosen  so  that  the  resulting  initial  ^-,,  after 

the  hurricane,  f  ,    was  subtracted  out,  was  as  smooth  as  possible 

in  the  vicinity  of  [x  (o),  y  (O)] .   The  calculated  hurricane 

trajectories  are  shown  in  Fig.  1.   The  trajectories  with  non-zero 

7,  viz.,  Tj/     interacting  nonlinear ly  with  f^,    always  move  to  the 
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left  of  the  trajectories  with  7=0  (non- interacting,  as  previously 
assumed  by  meteorologists).   However,  comparing  with  the  actual 
hurricane  trajectory,  there  are  certain  deficiencies  of  the  model 
(with  7  ?^0):   1)  the  model  is  slower  than  the  actual  motion,  2)  for 
the  initial  data  for  l4  August,  1955,  the  vortex,  laoving  initially 
in  a  northwesterly  direction  does  not  change  its  direction  to  the 
northeast  as  it  should.   The  run  with  7=0  does  make  the  turn 
(called  "recurvature"  by  meteorologists)  but  not  rapidly  enough. 
The  problem  of  predicting  recurvature  is  still  unsolved.   The 
vortex  trajectories  simulate  the  actual  trajectories  best  after 
recurvature;  and  the  time  development  of  ^^   is  also  good  for  these 
runs. 

(9) 
Interacting  Motion  of  Rectilinear  Qeostrophic  Bessel  Vortices'^'' 

In  (8)  and  in  the  problem  of  calculating  hurricane  trajec- 
tories, we  considered  only  a  single  geostrophic  vortex.   Now  we 
study  the  interacting  motion  of  N  (or  (N+l))  vortices.   The  equa- 
tions of  motion  can  be  derived'^'  from  (7).   For  N  vortices  at 
positions  (x.,y. ),  the  stream  function  satisfies 

P      N     ^  ^ 

(11.1)  (^  ''<)i'   =  JII  5(|r-r^|) 

and  the  solution  of  (11. l)  is 

(11.2)  '^  =  -^g^^i^of^l^-^il^ 

where  |r-r.  |  =  [(x-x.)^+  (y-y^)^]"^'^^.   If  the  only  elements  in  the 
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entire  flow  field  are  the  N  vortices,  the  velocity  of  the  k-th 
vortex  is  obtained  by  differentiating 

^11-5)         ^K)  =-^|^^iKo(-l^-^ol) 

i*l 

which  is  the  non- singular  part  of  (11,2)  at  (Xj^>yL-)>  yielding 

1  =  1        ^   ^ 

and 

1=1    ^  ^ 

Equations  (11.4)  and  (11,5)  are  a  coupled  system  of  2N  first 
order,  ordinary  differential  equations  for  [x.  (t  ),yj^(T )] ,  given 
the  initial  positions,  [x,  (o),y,  (o)]  .   These  equations  can  be 
expressed  in  a  more  symmetric  form, 

where 

(11.7)         ""^^^    ^iVc('<^lV^jf)  • 

Kirchhoff^  '    first  noticed  this  form  for  k   =   0.      Equations  (11.6) 
and  (11.7)  are  reminiscent  of  Hsunilton's  equations  of  motion  of 
interacting  particles;  however,  significant  differences  are  that 
1)  these  equations  are  in  the  two  dimensional  configuration  space. 
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{y^w>Yw)>    whereas  Hamilton's  equations  are  in  the  full  configura- 
tion-momentiim  phase  space,  and  2)  W  can  be  positive  or  negative 
depending  on  the  sign  of  the  r^'s.   The  integral  invariants  of 
Equations  (11.6)  and  (11.7)  are 

(11.8)  W  =  const.  ($  0)  , 

N  N 

(11.9)  m  r,^x  =  const.,   YH  ^k^k  =  const.  , 
k=l  ^  ^  k=l  ^  ^ 

N     p 

(11.10)  y"  7.  r  =  const.  , 

fer   ^  ^ 

(11-11)   g  VA  ^  7i7  ^  ^i^k' V^k'  %['c|r,-rj]  =  const. 

i,k=l 

where  (x^,y^)  =r^Jcos  6^^,  sin  6^^),  and  (*)  means  differentiation 
with  respect  to  time.   Equations  (11.9)  are  related  to  the  centroid 
of  the  vortices;  and  (11. lO)  and  (11. H)  are  reminiscent  of  momen- 
tum invariants.   For  k =  0,  (11.11)  reduces  to  the  classical  result: 

N       g.       N 

(11-12)  ^y^w\=^  Vk  • 

i,k=l 

We  study  the  linear  and  nonlinear  stability  of  the  inter- 
acting, kinematic  motion  of  (N+l)  gees trophic  vortices. 
Initially  the  vortices  are  positioned  near  a  uniformly  rotating 
equilibrium  configuration  which  consists  of  N  vortices  of  equal 
strength,  y,  equally  spaced  on  a  circle  of  radius,  a,  and  one 
vortex  of  strength,  y  ,  at  the  center  of  the  circle.  All  (N+l) 
vortices  have  the  seme  horizontal  length  scale  kT   ,  where  the 
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limiting  case,  fc  =  0,  corresponds  to  potential  (logarithmic) 
vortices.   The  ranges  of  the  parameters  N,  (7  ^7),  and  (ica.)    are: 
N  >  2,  -00  <    {■y   /y )    "^   00 ,    and  0  <_   (fca)  <  00  .   We  normalize  the 
continuous  parauneters  by  setting  7  =  1  and  a  =  1  to  simplify  the 
notation.   For  the  equilibi'ium  motion,  r,  =  1  (r  =  O)  and 
0   =  Ot  +  -i^  (e   =  o,  const.  )  where 

1    ^ 
(12)  ^  =  ^0"%^"^  "^  2  ?^ '^ki^i^'^ki^ 

i,k=l 

rr.fn                      M 1/2    -,       /      N    27r(k-i) 
where  a^^   =  /c[2(l-cos  w^^^)]  ''   and  w^^^  =  (Wj^-w^j  =  ^-^ '- <. 

For  the  linear  stability  analysis,  a  mixed  coordinate  representa- 
tion of  the  equations  of  motion  is  appropriate,  viz.,  polar 
coordinates  for  the  circle  vortices  and  Cartesian  coordinates  for 
the  center  vortex.   The  linearized  equations  of  motion  have  con- 
stant coefficients  in  a  frame  of  reference  rotating  with  angular 
velocity,  O.   Then  the  solution  ( r,   ,6^   ,  x *:   ,y);  ')    can  be  sought 

K      K      O      O 

A.t 
that  has  time  dependence  in  exponential  form  e    ;  and  if  an  eigen- 
value A.  has  a  real  part  that  is  positive,  the  motion  becomes 
exponentially  unstable  with  respect  to  the  equilibrium  motion.   The 
results  are  summarized  by  the  neutral  stability  curves  in  the 
(7  , /c )  parameter  space;  in  Fig.  2  a  circle  vortex  is  stable  on  the 
unhatched  (right)  side  of  each  curve;  in  Fig.  3  the  center  vortex 
is  stable  on  the  unhatched  (left)  side.   In  particular.  Table  I 
gives  the  7  -range  of  linear  stability  for  /c  =  0.   But,  a  complica- 
tion arises  because  of  the  presence  of  a  double  zero  eigenvalue 
for  all  values  of  the  parameters  N,  7   and  /c.   This  difficulty  is 
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partially  overcome  by  the  use  of  two  linear  integral  invariants; 
l)  by  linearizing  (11. lO)  we  get 


"    (1)  _ 


(«.l)  ^r^. 


=   const.  =  0  , 


and  2)  upon  linearizing  (11.11)  we  get 

(13.2)  "TZ   ^k   =  const.  ^  0 

x^rhich  cannot  be  set  equal  to  zero  usually.   Since  these  first 
order  (linear)  invariants  do  not  include  the  center  vortex 


position  (x^-^Sy^'*''),  which  only  appears  to  second  (and  higher) 
order,  the  double  zero  eigenvalues  can  be  strictly  nullified  only 


for  7   =  0.   For  7  7^  0,  the  double  zero  eigenvalues  allov;  solu- 
tions which  increase  linearly  with  time,  implying  that  the  motion 
of  the  center  vortex  is  nonlinear.   Numerical  integrations  of  the 
nonlinear  equations  of  motion  viere  carried  out.   The  numerical 
computations  show  that  nonlinear  effects  are  important  in  a 
relatively  small  transition  region  embracing  both  the  stable  and 
unstable  sides  of  the  neutral  stability  curves  for  the  circle 
vortices  (Fig.  2),  even  for  small  perturbations  of  the  initial 
position;  othervrise  Fig.  2  represents  valid  results  concerning 
the  stability  of  the  circle  vortices.   However,  the  stability  of 
the  center  vortex  is  not  well-represented  by  the  neutral  stability 
curves  sho\m  in  Fig.  3 i    and  the  numerical  computations  are 
essential  for  the  study  of  the  stability  of  the  center  vortex. 
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(13)  \    114)  115]  (:3)  liri  1161  (19)  120)  121)  1221  1231 

Fig,   1.     Gcostrophic  point  vortex  trajectories  for  hurricane  "Betsy"   initial  data. 
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Fig.  2  Neutral  stability  curves  for  circle  vortices  (stable  region      Fig. 3  .  Neutral  stability  curves  for  center  vortex  (stable  region 
on  right  side  of  each  curve) .  on  left  side  of  each  curve) . 

Table  I.  Range  of  exponential  stability  for 
free  center  vortex  case,  k  =  0. 


N 


Lower  stability        Upper  stability 
limit*  limit'' 


2 

7o'  =  —  " 

7o"=-1.2S 

3 

—  CO 

1 

4 

-0.5 

2.25 

S 

-O.S 

4 

6 

-0.25 

6.25 

7 

0 

9 

8 

0.5 

12.25 

9 

1 

16 

10 

1.75 

20.25 

11 

2.5 

25 

12 

3.5 

30.25 

13 

4.5 

36 

14 

5.75 

42.25 

IS 

7 

49 

*  A  circle  vortex  is  unstable  for  70  less  than  the  lower  stability  limit. 
^  The  center  vortex  is  unstable  for  70  greater  than  the  upper  stability 
limit.       .    ,  , 
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Lecture  9 

The  Upper  Boundary  Condition  for  Gravity  Waves  in  the  Atmosphere 

Michael  Yanowitch 

The  linearized  theory  of  waves  in  an  inviscid  isothermal 
fluid  in  a  half  space  serves  as  a  simple  model  for  the  study  of 
atmospheric  gravity  waves.   Two  boundary  conditions  are  appropriate 
for  this  model:   a  boundary  condition  on  the  vertical  velocity  at 
the  ground,  and  an  "upper  boundary  condition",  which  depends  on 
the  form  of  the  solution.   If  the  solution  is  nonoscillatory  in 
the  vertical  (z)  coordinate,  one  imposes  the  condition  that  the 
kinetic  energy  in  an  infinite  vertical  column  of  fluid  should  be 
finite.   If  the  solution  is  wavelike  in  the  x-coordinate,  a  radia- 
tion condition  is  customarily  imposed,  i.e.  it  is  assumed  that 
the  energy  propagates  upward  to  infinity  without  reflection. 
Unfortunately,  the  validity  of  this  model  is  open  to  question 
since  the  velocities  increase  exponentially  with  altitude  and  thus 
violate  the  assumptions  underlying  the  basic  linearization.   To 
overcome  this  deficiency  we  will  study  a  viscous  fluid  model.   It 
will  be  seen  that  the  solutions  to  this  problem  are  uniformly 
small,  consistent  with  the  small  amplitude  assumption,  and  that 
the  radiation  condition  is,  in  general,  not  correct  since  waves 
may  be  reflected  downward.   For  our  purpose  it  suffices  to  consider 
an  incompressible,  but  horizontally  stratified,  fluid. 

The  undisturbed  incompressible  fluid  occupying  the  half-space 
z  >  0  is  assumed  to  have  a  constant  dynamic  viscosity  e,  and  a 
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stratified  density  p(z),  a  specified  function.   The  wave  motion  has 
a  velocity  u  in  the  horizontal  direction  x  and  a  velocity  w  in  the 
vertical  direction  z,  with  a  density  perturbstion  p  and  a  pressure 
perturbation  p.   The  linearized  differential  equations  governing 
the  motion  are         .  ,     i    ■,,. 

^  u  +  ^  w  =  0 
ox      cz 


S  ~ 


-TT  P  +  5^  W  =  0 

Zz   ^        dz 


P  3t 

p  Zt 


u  +  -^  p  =  e 
ox  ^ 


,\2  ^2  \ 

\  ^.x^     az  V 


w  + 


^  P  +  gp  =  e 


where  g  is  the  acceleration  of  gravity.   The  stream  function  ^,  so 

defined  that  u  =  -r^  and  w  =  -  -y—,    satisfies 

oz  ox 


S' 


dt 


n2       n     V 

\~2  oz  ^    £z 

ox 


-  ai£ 


^^^ 


T  =  e 


upon  elimination  of  the  density  and  pressure  perturbations. 

For  a  wave  propagating  in  the  x-direction,  we  look  for  a 
solution  such  that  T  =  ^{z)e^   x-icu  ^  ^^^   complex-valued  amplitude 
Ci(z)  then  satisfies  ,   ,  L' 


-J—  p  -r-  <f>  +  K 

dz  ^   dz 


p-^H 


0  =  0. 


0) 


The  no-slip  condition  at  the  lower  boundary  requires  that  u  =  0 
and  w  =  const,  e"""  ^"^'^  at  z  =  0,  i.e. 
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<I)  I    =  const.  ,   -r— 
lz=0  clz 


0  . 


Since  the  density  in  the  atmosphere  decreases  very  rapidly,  by  a 
factor  of  10  in  the  lowest  100  km.  (while  the  viscosity  varies 
only  slightly),  we  shall  assume  that  the  density  is  stratified 
exponentially,  viz.. 


p(z)  .  p^e-^/H 


where  H  is  a  scale  height  for  the  stratification.  Accordingly,  0 
satisfies 


dz 


-  k'"  I  $  +ia)p  e 


•z/H 


dz  V  CO  H 


-  1  U 


=  0  . 


Instead  of  the  radiation  condition,  we  shall  impose  a 
"dissipation  condition":   the  average  rate  of  energy  dissipation 
in  an  infinite  vertical  column  of  fluid  (O  <  z  <  oo  )  of  unit 
cross-section  area  must  be  finite.   Namely, 


00  00  00^  2  \ 

0  0  0     ^^^   / 


dz  <  00  , 


as  the  dissipation  function  depends  on  the  square  of  the  space 
derivatives  of  velocities. 

At  very  high  altitudes  (z  -*■  oo  )  the  viscous  term  dominates, 
Hence  in  this  "viscous  region" 

'd^     2Y 

dz'^    y 
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thus,  $  is  a  linear  combination  of  e    and  ze    in  order  to 
satisfy  the  dissipation  condition.   Separated  from  the  viscous 
region  by  a  "transition  region"  (where  e"  e'^^  =  0(l))  is  an 
"inviscid  region",  where  the  viscous  effect  is  negligible,  hence 


thus 


with 


e.  ^   Ae^/2H-iP^  +  Be^/2«-^i^2 


P  =  [k2(gH/a)2  -  H^)  -  ^]^/^ 


The  inviscid  region  is  connected  to  the  lower  boundary  z  =  0  by  a 
"boundary  layer"  which  reduces  -^   to  zero  as  z  -*  0.   The  problem 
is,  therefore,  to  find  a  solution  which  is  valid  uniformly  in  the 
inviscid,  transition,  and  viscous  regions.   (The  boundary  layer 
can  always  be  fitted  in  later. ) 

To  obtain  the  solution  it  is  advantageous  to  transform  to  the 
independent  variable:  i  =   (-i/e)e"^^  .   The  transformed  equation 
(in  dimensionless  quantities)  is 


where  r  is  a  constajit.   The  only  singular  points  of  the  differen- 
tial equation  are  |  =  0  and  ^  «  oo  ,  with  ^  ■=  0  being  a  regular 
singular  point,  ^  =  oo  an  irregular  one.   It  is  possible  to  obtain 
the  unique  solution  satisfying  the  dissipation  condition  in  the 
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form  of  a  Frobenius  series  about  ^  =  0.   The  coefficients  of  the 
series  satisfy  a  two  term  recursion  formula,  which  enables  one  to 
convert  the  series  to  a  complex  integral  by  the  Cauchy  summation 
method.   The  behavior  in  the  inviscid  region  as  e  ->•  0  corresponds 
to  the  behavior  for  large  ^,  with  arg  ^  =  -Iv/k,    which  can  be 
obtained  by  evaluating  the  residues  at  the  poles  of  the  integrand. 
The  solution  does  not  tend  to  a  limit  as  e  -*■  0  in  the  fixed 
X,  z  coordinate  system,  but  does  so  in  a  coordinate  system  which 
moves  upward  as  e  ~*  0  in  such  a  way  that  e~  e"^       =   const.   The 
magnitude  of  the  reflection  coefficient,  fB/A|,  however,  tends  to 
the  limit  exp  {~2k   H/x),  where  ^  is  the  vertical  wavelength.   Thus 
there  is  substantial  reflection  when  the  vertical  wavelength  is 
large  and  the  solution  approaches  the  form  of  a  standing  wave  as 
X  -«•  oo .   The  radiation  condition  becomes  more  and  more  accurate 
as  >v  -f  0. 
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Lecture  10 

Numerical  Simulation  of  the  Life  Cycle 
of  Tropical  Cyclones 

Katsuyuki  Goyama 

The  tropical  cyclone  is  a  solitary  creature  of  the  tropical 
oceans  accompanied  by  violent  rotating  winds  and  torrential  rain. 
Observational  studies  and  diagnostic  analyses  leave  little  doubt 
that  the  energy  required  for  driving  the  vortex  comes  from  the 
latent  heat  of  condensation  released  by  tall  convective  clouds 
around  the  center,  and  that  the  frictionally  induced  inflow  in  the 
vortex  plays  a  major  role  in  supporting  the  continued  activity  of 
convective  clouds.   This  dual  character  with  respect  to  important 
scales  of  motion  poses  a  great  difficulty  in  investigating  the 
dynamics  of  tropical  cyclones  as  time-dependent  phenomena.   How- 
ever, in  order  to  understand  the  large-scale  aspects  of  tropical 
cyclones,  one  may  formulate  the  role  of  convective  clouds  in  terms 
of  cyclone-scale  variables  with  only  implicit  consideration  of  the 
dynamics  of  individual  clouds.  The  present  study  is  an  attempt  to 
understand  the  basic  mechanism  of  tropical  cyclones  by  constructing 
a  numerical- dynamical  model  on  such  a  basis. 

The  model  assumes  that  the  large-scale  hydrodynamlcal  aspects 
of  a  tropical  cyclone  may  be  represented  by  an  axisyrametric,  quasi- 
balanced  vortex  in  a  stably  stratified  incompressible  fluid,  while 
the  effects  of  moist  convection  may  be  formulated  through  the  first 
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law  of  thermodynamics  applied  to  an  implicit  model  of  penetrative 
convective  clouds.   The  air-sea  exchange  of  angular  momentum  as 
well  as  latent  and  sensible  heat  is  explicitly  calculated  in  the 
model  with  the  use  of  conventional  approximations. 

Results  of  numerical  integration  show  that  the  model  is 
capable  of  simulating  the  typical  life  cycle  of  tropical  cyclones, 
including  the  mature  hurricane  stage,  with  a  remarkable  degree  of 
reality.   The  response  of  the  model  cyclone  to  changes  in  such 
parameters  as  the  sea  surface  temperature,  the  coefficient  of 
air-sea  energy  exchange,  and  the  initial  conditions  is  tested  in 
a  number  of  numerical  experiments  to  show  quite  plausible  results, 
A  detailed  diagnosis  of  the  energy  budget  of  the  simulated 
tropical  cyclone  is  also  carried  out.   The  rate  of  total  rainfall, 
the  production  and  dissipation  of  kinetic  energy,  and  other 
energetic  characteristics  of  the  computed  cyclone  compare  very 
well  with  available  estimates  for  observed  tropical  cyclones. 

Because  of  the  restrictive  assumption  of  axisymmetry  and 
other  weak  approximations,  the  model  is  not  realistic  enough  to 
predict  behavior  of  individual  tropical  cyclones  in  nature.   The 
limitation  of  the  present  model  in  this  regard  is  also  discussed. 

Figure  1  shows  the  characteristic  scales  (lengths,  times, 
and  speeds)  for  various  meteorological  phenomena. 
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Characteristic    Time    Scale    T   in   Seconds 
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Lecture  11 

A  Mountain  Building  Model 
(A  presentation  of  the  Ph.D.  Thesis  of  Jay  Wolkowisky) 

Chester  B.  Sensenig 

The  mountain  building  model  treated  in  the  Ph.D.  thesis 
of  Jay  Wolkowisky  [1]  was  proposed  by  J.  J.  Stoker.   It  is 
based  on  the  observation  that  in  a  number  of  instances  relatively 
high  mountain  ranges  and/or  deep  oceajn  troughs  occur  near  the  edge 
of  a  continent,  and  the  continent  is  relatively  flat  except  for 
these  variations  in  elevation  near  its  edge.   Since  the  earth's 
crust  over  a  continent  is  much  thicker  than  over  most  of  the  ocean 
floor  but  is  thin  relative  to  the  radius  of  the  earth,  the 
continent  can  be  thought  of  as  a  thin  shell  with  its  edge  at  or 
near  the  ocean  shore.   Stoker  then  posed  the  question, 
"Can  the  relatively  flat  interior  and  rippled  edge  of  the  continent 
be  explained  as  the  lowest  buckled  mode  of  a  thin  elastic  plate 
resting  on  an  elastic  foundation?"   For  simplicity,  a  thin  circular 
plate  resting  on  an  elastic  foundation  is  considered,  and  a 
compressive  edge  force  is  applied  with  the  edge  clamped,  or 
pinned  (zero  bending  moment).  Values  of  foundation  stiffness, 
plate  thickness,  and  edge  thrust  are  sought  for  which  the  lowest 
buckled  mode  is  relatively  flat  in  the  interior  with  relatively 
large  ripples  near  the  edge. 
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The  buckling  problem  is  treated  by  using  the  plate  theory 
of  V.  Karman  and  Fbppl  [2]  with  a  term  added  which  corresponds 
to  the  elastic  foundation.   The  differential  equations  are 

■y"h  EAw  =  (l)W   -2(1)  w   +cl)W   -trW 


Vy  XX   ^xy  xy  ^xx  yy   h 


1   A^i        2 

^Acp  =  w   -w  w 
^   ^    xy   XX  yy 


2    2    2 
X  and  y  are  the  independent  variables  with  x  +  y  _<  R 

(R  =  the  radius  of  the  plate), 

vj  =  w(x,  y)  is  the  vertical  deflection, 

(j)  =  (f)(x,  y)  is  the  Airy  stress  function, 

T   =(j),T   =-<{),  T   =  <i>        are  the  longitudinal  stresses 
XX   ^yy'   xy    ^xy'   yy   ^xx  ^ 

averaged  through  the  thickness  of  the  plate, 

2       1 

■y     =   __   (v  =  Poisson's  ratio), 

12(l-v'') 
E  =  Young's  modulus, 

h  =  the  plate  thickness, 

k  =  the  foundation  stiffness  (kw  =  restoring  force  per  unit 

area), 
>2    ^2 

^x^   ^y^ 
and  the  subscripts  x  and  y  indicate  partial  derivatives  except 

in  the  case  of  the  longitudinal  stresses. 
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In  the  books  of  Jeffry  (The  Earth,  3rd  Ed. ),  and  Helskanen 

(The  Earth  and  Its  Gravity  Field),  it  is  stated  that  elastic 

instability  of  the  earth's  crust  when  under  compressive  stress 

never  can  occur  because  the  stress  required  would  be  much  greater 

than  the  crushing  strength  of  the  material  of  the  crust.   This 

follows,  however,  only  because  these  writers  assume  buckled  forms 

for  the  crust  of  the  earth  that  have  very  small  wavelengths,  of 

the  order  of  roughly  three  times  the  thickness  of  the  crust.   They 

use  the  Euler  theory  of  buckling  of  columns  in  their  discussion, 

but  this  theory  applies  only  for  long  slender  columns  (having 

lengths  of  the  order  of  at  least  50  times  the  thickness,  say)  and 

gives  far  too  large  values  for  thick  columns «   In  addition  there 

seems  no  reason  not  to  investigate  the  elastic  stability  of  the 

crust  under  compression  when  it  is  assumed  to  be  a  plate  of  the 

size  of  a  continent.   There  is,  for  exeonple,  a  geosyncline  in  a 

large  portion  of  the  interior  of  continental  North  America,  with 

a  wavelength  of  something  like  I50O  miles.   It  is  assumed  here 

that  such  a  plate  buckles  under  compression  into  a  cylinder  with 

generators  in  the  north-south  direction,  so  that  (l)  applies  when 

w  and  (\)   depend  only  upon  x.   It  is  assumed  further  —  as  was  done 

by  the  two  writers  mentioned  above  —  that  the  ends  of  plate  on  the 

east  and  west  sides  are  free  of  bending  moments,  which  implies  that 
2 

— 7y  =  0  there,  and  that  the  effect  of  the  support  of  the  earth 

dx 

■under  the  crust  can  be  ignored. 

Only  the  onset  of  buckling  is  to  be  determined,  and  that  can 

be  done  by  linearizing  (1)  (as  will  be  explained  in  more  detail 

later),  and  setting  i>        =   -t      ,    i)        =   i)        =Oin  which  t    is  the 
"  ^   ^yy     xx'  ^xy    ^yy  xx 
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compressive  stress  assumed  to  exist  in  the  x-direction  (east-west). 
The  weight  of  the  earth  is  neglected.   The  result  is  the  following 
problem.   The  differential  equation 

2     4        -2 
h  E    dw^_^    dw^    0   <  X   <   £    , 

12(l-v^)  d?     ^^  d?  '      "   ~   ' 

is  to  be  integrated  under  the  boundary  conditions  w  =  0,  and 
.2 

— ^  =  0  at  X  =  0  and  x  =   H,  o      This  is  a  linear  eigenvalue  problem 
dx 
for  the  determination  of  the  buckling  stress  t   ,  i.e.  the  smallest 

value  of  T    for  which  a  nonvanishing  solution  of  the  mathematical 

XX 

problem  exists.   The  solution  is  well  known  and  is  given  by 

TT      h  E 


"^xx   -J2   i2(l-v2) 

If  V  =  0.3,  i    =  1500  miles,  E  =  12x10   lbs/in  ,  and  h  =  21  miles 
are  taken,  the  value  of  r   is  found  to  be 

2/1  n2  12x  10^ 
^xx  "^  ^  ^75^   12(0.9) 

~  2100  lbs/in^ 


which  is  negligible.   For  £   =   1000  miles,  for  which  the  theory 

'I  ~  50),  T 
■h   -^  ^'      XX 


would  still  be  valid  (^  ~  50),  t   -^   4800  lbs/in  holds.   Even  for 


I   =  500  miles,  the  buckling  stress  is  only  I90OO  lbs/in  .   These 
values  are,  of  course,  far  under  the  crushing  stress  for  materials 
of  the  type  making  up  the  crust  of  the  earth.   Thus  it  would  seem 
not  reasonable  to  reject  buckling  of  the  earth's  crust  due  to 
elastic  instability  as  a  possible  effect  in  interesting  special 
cases. 
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For  simplicity,  radially  symmetric  deformations  of  the  plate 
are  considered.   Thus  the  partial  differential  equations  (l) 
become  ordinary  differential  equations. 

Thesa  equations  (with  k  =  O)  were  used  by  Friedrichs  and 

Stoker  [J]  to  study  the  lowest  buckled  mode  of  a  circular  plate 

due  to  a  compressive  edge  force  with  pinned  edge.   It  was 

discovered  that  the  solution  depends  essentially  on  N  =  -^  where 
_  p 

p  is  the  pressure  applied  at  the  edge  and  p°  is  the  lowest 

pressure  at  vrhich  the  plate  buckles.   The  solution  was  calculated 

for  1  <  N  <  2.5  by  using  a  perturbation  method,  for  2.5  ^  N  ^  15 

by  using  power  series,  and  for  N  -*  00  by  using  asymptotic  methods o 

Higher  buckled  modes  were  studied  by  others  including 
H.  Keller,  J.  Keller,  and  Reiss  [4]  and  Conn  [5] . 

From  these  treatments  it  is  known  that  the  plate  is  rela- 
tively flat  in  the  interior  with  large  variations  in  w  near  the 
edge  if  the  compressive  edge  force  is  sufficiently  large.   The 
lowest  mode  of  buckling  will  look  like  a  plain  with  a  cliff  at 
the  edge,  but  the  higher  modes  of  buckling  will  have  hills  and 
valleys  near  the  edge.   However,  the  higher  modes  of  buckling 
certainly  will  not  occur  in  nature.   The  question  being  raised 
is  whether  the  lowest  buckled  mode  will  be  relatively  flat  in  the 
interior  and  have  hills  and  valleys  near  the  edge  for  some 
parameter  values  when  the  plate  is  on  an  elastic  foundation. 

The  following  quantities  are  introduced: 
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R  =  radius  of  the  plate. 


r*  =ix2-Hy2, 


p  =  edge  thrust, 

/   * \  Fe  R  dw 

p  =  -  2_  (mean  radial  stress), 

pr*  dr 

* 
r 

W  =  -  -  -^^   /    tw(t)dt, 
4   p  Rr*2  ^Q 

7  h  E 


7  h  E 

1  d  ^3  d  s 
^  -  p  d7  ^"^  d?^' 

The  differential  equations  then  become 

^  Gq  +  a^pq  =  K^W 
(2)  J   Gp  =1  q2 

1^  GW  +q  =  0  . 

To  obtain  (2)  the  equations  (l)  are  integrated  and  integra- 
tion constaxits  a^re  eliminated  by  using  regularity  and  symmetry  of 
stress  and  displacement  at  the  center  of  the  plate. 

The  boundary  conditions  are: 
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(3) 


f   q'(0)  =  p'(0)  =  W'(0)  =  0  from  symmetry, 
p(l)  =  1  (edge  pressure  is  p), 
W'  (1)  +2W(1)  =  0  (w  =  0  at  edge), 
q'(l)+  (l+v)g(l)  =  0  for  pinned  edge, 
q(l)  =  0  for  clamped  edge. 


It  is  seen  that  the  problem  depends  essentially  on  the 
parameters  a  and  K. 

A  perturbation  method  is  used  to  determine  the  critical 
value  of  the  parameters.   Since  an  exact  solution  to  (2)  and  (3) 
is  p  =  1,  q  =  W  s  0,  we  assume  an  expansion 

2 

q  =  eQi  +  e  qo+  •  •  • 


P  =  1+  ePT+  e  Po  + 


2-  . 


W  =  eW-L+  e  W2+  .. 


where  e  is  a  perturbation  parameter. 
From  (2)  and  (3)  we  obtain 


(^) 


and 


Gq.j^+  a^q^  =  K  W^  , 


Gp^  =  0 


GWtl+  q-j^  =  0 
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^ci[{0)   =  Pi(0)  =  W[(0)  =  0  , 
Pl(l)  =  0  , 


(5)   /   W|(1)+2W^(1)  =  0  , 

q^(l)  +  (l+v  )q-,^(l)=  0  for  a  simply  supported  edge. 


q-,(l)  =  0  for  a  clamped  edge. 


These  imply  p-,  =   0,  and  leave  two  differential  equations  for 

2 
q,  euid  W-,.   Special  solutions  can  be  found  such  that  q-j^  =  A  Wj^ 

where  A  is  a  constant.   From  (4) 


so  that 


GWj^  :-  >v^W^  =  0  , 


^2    2   K"" 


=  0  , 


The  solutions  A  are 


H  = 


^-^ 


1  ry 


hK^  , 


(6) 


1^-1/ 


?Cte5 


where  we  are  nov/  limiting  ourselves  to  pa-ametor  values  such  that 


(7) 
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so  that  >^n  and  "Kp   are  real. 

The  general  solution  for  q^  and  W^  is 

W-L  =  i[AJ3^(X3^r)  +  BY^(X^r)  +  CJ3^(>y2r)  +  DY3^(7.2^)]  , 

^1  "  -^tA^iJiCT^ir) +B?y^Y3^(^J^r)  +  C^|J^(7^2r)+ DTvgY^Cx^r)]  , 

where  J,  and  Y,  are  the  usual  Bessel  functions  and  A,  B,  C,  D  are 
arbitrary  constajits. 

If  we  substitute  W,  and  q,  into  the  boundary  conditions  (5), 
we  obtain  4  homogeneous  linear  equations  in  A,  B,  C,  D.   In  order 
that  there  be  a  non-zero  solution^  the  determinant  of  the 
coefficients  must  vanish.   This  gives 

\   T-Orrr  -  T:7  ^  ^2  T^JTf   -  HV  f°^  ^  Pi^ed  edge  , 


'o^'^l'   "^    "  *'o^"2 


(8) 


^1  J  (\ — r  ~  "^2  TTx — r  ^°^  ^  clamped  edge  . 


Equations  (8)  determine  the  relationship  between  a  and  K  at 
the  onset  of  buckling.   For  each  K  there  will  be  a  sequence  a  (k) 
of  values  of  a  such  that  (8)  is  satisfied.  We  let  A,  (K)  and 
^2n^^^  be  the  values  of  7,^   and  >,„  which  correspond  to  o,f.M   and  K. 
We  also  let  q,   and  W,   be  the  functions  q,  and  W,  corresponding 
to  a^(K),  K,  and  we  let  -w-^n   ^®  ^^^  egrresponding  perturbed 
vertical  deflection.   Then 
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(9)  ;  Wi^=^[^2nJo<^2n>Jl(^ln^)-^ln'^o^Hn^Jl^^2n^^^  ' 


V  ^ln  =  VJo(^2n)Jo(^ln^^-'^o^Hn)Jo('^2n^)l  > 


where  a  and  b  are  arbitrary  constants, 
n      n 

From  (8)  and  (6)  it  follows  that  as  K  -♦  0,  a^(K)  goes  to  the 
n^'^  critical  value  for  the  plate  without  an  elastic  foundation. 
Hence,  although  in  (7)  we  restricted  the  range  in  which  we  looked 
for  critical  values  of  a,  we  still  have  a  critical  value  of  a 
which  corresponds  to  each  of  those  when  the  elastic  foundation  is 
missir\g.   Since  the  elastic  foundation  is  expected  to  have  a 
stiffening  effect  on  the  plate,  we  believe  we  have  all  the 
critical  values  of  a  despite  the  restriction  (7). 

From  (9c) 

^11  =  ^lfJo(^21^Jo(^ll^)-  Jo^hl^Jo^'^21^)^  • 

Since  7^,,  -♦  co  as  K  -►  oo  ,  it  is  seen  that  w,,  can  be  made  to 
have  as  many  ripples  as  one  pleases  by  picking  K  appropriately, 
i.e.  the  lowest  mode  of  buckling  can  be  made  to  have  as  many 
ripples  as  one  pleases  by  picking  K  appropriately  and  then  picking 
the  edge  thrust  slightly  larger  than  a-,(K).  ,.  ■ 

Numerical  results  have  been  obtained  for  the  clamped  edge 
case.   Figures  1  and  2  show  A,  (K)  and  a  (K).   Figures  3  through  9 
show  the  dimensionless  vertical  deflection  for  a  variety  of  values 
of  a  and  K  and  for  n  =  1,2, 3. 
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By  studying  Pigs.  3-9  one  concludes  that  if  K  is  held  fixed 
and  a  is  increased,  the  zeros  of  w  disappear  or  move  near  to  the 
edge  of  the  plate  while  the  interior  of  the  plate  approaches  a 
flat  configuration.   In  particular,  all  zeros  of  the  lowest  buckled 
mode  disappear  as  a  -*  oo  with  K  fixed.   Consequently,  this  pro- 
cedure does  not  lead  us  to  a  solution  of  the  type  sought,  i.e. 
a  lowest  buckled  mode  with  a  relatively  flat  interior  and.   ripples 
near  the  edge. 

The  next  step  is  to  see  what  happens  if  we  let  a  -»■  oo  and 

K  -♦  00  together.  Since  the  critical  values  of  a  are  all  in  the 

2 
range  a  >  2K  and  since  we  want  a  to  go  to  oo  in  such  a  way  that  a 

is  slightly  greater  than  a,(K)  but  not  too  much  greater,  we  con- 

p  p  oh. 

sider  values  of  a  >2K  or  0  < K/a  < l/2.   Letting  k  =  K^/a  ,  this 

range  becomes  0  <  k  <  1/4.   It  is  natural  therefore  to  try  letting 

K  -►  00  with  k  a  constant  in  the  interval  (O,  ^).   Numerical  results 

were  obtained  using  various  values  of  k  for  the  first  and  second 

modes  of  buckling.   No  significant  qualitative  differences  were 

obtained  except  that  limiting  values  were  approached  more  rapidly 

as  K  ->•  00  if  k  was  small.   Figures  10  and  11  show  the  dimensionless 

radial  bending  stress  and  vertical  deflection  for  large  values  of 

a  and  K  with  k  =  .01  .  The  solution  shown  in  Fig.  11  is  of  the 

type  being  sought.  It  is  very  nearly  flat  with  many  small  ripples 

in  the  interior  eind  has  two  relatively  large  hills  or  mountains 

near  the  edge. 
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The  thesis  also  contains  a  formulation  and  discussion  of  the 
boundary  layer  problem  which  results  if  a  -*  oo  with  K  fixed  or  if 
a  -*  00  and  K  -♦  oo  with  k  fixed.   It  is  observed  that  in  the  former 
case  the  boundary  layer  problem  is  exactly  the  same  as  that  treated 
by  Friedrichs  and  Stoker.   It  is  concluded  that,  as  a  increases 
with  K  held  fixed,  the  plate  behaves  like  a  plate  without  elastic 
foundation.  The  boundary  layer  problem  corresponding  to  the  second 
case  could  not  be  treated  by  the  methods  used  by  Friedrichs  and 
Stoker,  and  it  remains  unsolved. 
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//4?F 
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dimenSionless  deflection 
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w//^R 
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DIMENSONLE'rS  DEFLECTION 
TYPE^O.2  SOLUTION,  K  =  I18,    03=  18.4 


I'i'  lire  7 
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0  =  2^26 
c=52.90 
a  =64.98 


a=52.9a 
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Lecture  12 

T\irbulence 
J.  B.  Keller 

At  low  velocities  the  fluid  motion  appears  as  a  laminar  flow. 
It  may  be  steady  or  it  may  change  slowly  if  the  external  condi- 
tions are  changed.  At  high  velocities  the  flow  becomes  turbulent. 
It  changes  rapidly  in  time  even  when  the  external  conditions  are 
constant.   For  each  flow  configuration,  there  is  a  critical  value 
of  the  Reynolds  number  R  =  UL/v,  (U  is  a  typical  velocity,  L  is  a 
typical  length  and  v  is  the  kinematic  viscosity)  such  that  the 
flow  is  turbulent  if  the  Reynolds  number  exceeds  the  critical 
value .  ' 

The  first  problem  of  the  theory  of  turbulence  is  to  account 

for  the  phase  transition  of  a  flow  from  the  laminar  state  to  the 

turbulent  state,  and  to  provide  a  procedure  for  determining  the 

transition  point,  i.e.  the  critical  Reynolds  number  R  .  The 

linear  stability  theory  offers  the  answer.  The  laminar  state  is 

stable  if  R  <  R„  and  unstable  if  R  >  R^.  Therefore  when  R  >  R„, 
c  c  c 

any  small  fluctuation  in  the  initial  velocity  of  the  flow  or  in 
the  velocity  of  a  boundary  of  the  fluid  will  grow  and  convert  the 
laminar  flow  into  a  turbulent  flow. 

A  second  problem  of  the  theory  of  turbulence  is  to  determine 
the  behavior  of  flows  with  R  slightly  larger  than  R..  These  flows 
can  be  viewed  as  laminar  flows  combined  with  a  few  unstable  modes 
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having  small  amplitudes.   Bifurcation  theory  shows  that  the  number 
of  steady  solutions  increases  with  the  Reynolds  number.  Thus, 
according  to  the  nonlinear  stability  theory  the  turbulent  flow 
may  be  regarded  as  approximating  one  solution  for  a  while,  then 
another,  etc. 

The  main  problem  of  turbulent  theory  is  to  describe  flows 
with  R  much  larger  than  R  .   The  formulation  of  a  theory  of  fully 
developed  turbulence  is  the  major  goal  of  turbulence  theory,  but 
that  goal  is  nowhere  near  being  achieved. 

Fully  developed  turbulence  occurs  in  various  fluid  flows  in 
geophysics  and  astrophysics.   The  lack  of  a  theory  of  such  flows 
is  a  major  barrier  to  progress  in  all  these  fields. 


Reference  ,i   ; 

J.  B.  Keller,  Survey  of  the  theory  of  turbulence.  Contemporary 

Physics,  vol.  I,  International  Atomic  Energy  Agency,  Vienna,  1969. 
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Lecture  13 

Progressive  Gravity  Waves  on  a  Sphere 
A.  S.  Peters 

A  remarkable  example  of  an  axisymmetric  disturbance  which 
propagates  on  the  surface  of  a  sphere  is  the  atmospheric  pressure 
wave  exhibited  during  the  eruption  V  of  the  volcano  Krakatoa  in 
1883.   In  that  event,  the  disturbance  originated  at  a  pole  point 
(Krakatoa).   It  propagated  over  the  surface  of  the  earth  almost 
axisymmetrically.   Then  the  pressure  wave  was  reflected  at  the 
antipodal  point  located  in  Venezuela,   Several  passages  of  the 
disturbance  were  recorded. 

Progressive  waves  on  a  sphere  can  be  studied  using  a  linear 
model  with  an  incompressible,  inviscid,  gravitating  fluid.   In 
spherical  coordinates  (r,  0,(j))  the  elevation  of  the  free  surface 
^(e,t)  satisfies  the  wave  equation  on  a  sphere 

,2 


1    S  ^.   o  S   ,     1     ^'^ 


where  c  =  /gh/a,  g  being  the  gravitational  acceleration,  h  being 
the  depth  of  the  shallow  fluid,  and  a  being  the  radius  of  the 
sphere.   The  initial  conditions  to  be  satisfied  are 

^1    =  0  ,    1^1    =  f(0) 
where  f(e)  is  the  initial  disturbance. 
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Formally  the  solution  can  be  obtained  as  an  eigenfunction 
expansion  in  Legendre  polynomials.   But  the  intrinsic  nature  of 
wave  propagation  is  not  revealed  by  an  infinite  series  of  func- 
tions involving  the  time.   However,  by  applying  the  Laplace  trans- 
form with  respect  to  the  time  the  solution  can  be  expressed  as 
the  sum  of  a  finite  number  of  terms.   In  this  way  a  clearer  picti.^-e 
of  the  wave  motion  can  be  obtained.   For  detailed  calculations 
see  IMM-NYU  Report  271,  June  I960,    VJaves  on  a  sphere,  by 
A.  S.  Peters. 


V "Eruption  of  Krakatoa  and  Subsequent  Phenomena".   Report  of  the 
Krakatoa  Committee  of  the  Royal  Society  of  London  (1888). 
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Lecture  l4 

Turbulence  Spectra  and  Vortex  Formation 
Alexandre  Joel  Chorin 

Introduction.  It  is  well  knovm  that  the  energy  spectnim 
E(k, t)  of  a  turbulent  flow  has,  for  large  wave  number  k,  a 
universal  form,  i.e.  a  form  common  to  a  variety  of  flows  whose 
larger  scale  features  have  little  in  common.  This  fact  is  often 
explained  by  the  universal  equilibrium  hypothesis,  well  summarized 
in  Batchelor  :   The  range  of  wave  numbers  k  which  contain  most  of 

the  energy  ("the  energy  containing  eddies")  can  be  regarded  as  a 

""?  1/2 
definite  group,  with  characteristic  velocity  u  =  (u  )  '   and 

Characteristic  length  2   »  k^^,  where  u  is  the  velocity  vector,  the 
bar  under  u  denotes  a  vector,  the  bar  above  (u  )  denotes  an 
appropriate  average,  and  k   is  a  typical  wave  number  in  the 
group.  The  characteristic  time  of  these  eddies  is  l/n,    and  the 
time  scale  of  their  decay  is  u/|3H|;  these  times  are  experimentally 
found  to  be  comparable,  and  thus  this  range  of  k's  has  no  feature 
resembling  an  equilibrium.  It  may  however  be  assumed,  under 
conditions  to  be  determined,  that  for  large  k  the  spectrum  has  a 
characteristic  time  small  in  comparison  with  the  scale  of  the 
over-all  decay,  and  thus  will  be  associated  with  degrees  of  free- 
dom in  approximate  statistical  equilibrium.   Let  u(k)  be  the 
Fourier  transform  of  u(x),  let  k   be  a  wave  number  typical  of  the 
equilibrium  range,  with  magnitude  k^  =  |k^q U  and  let  u   be  a 
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typical  amplitude  of  u(k),  k  =  |k|  in  the  equilibrium  range,  for 
example  u^^  -  t5(]Seq^l*   ^* 

The  characteristic  time  of  u^  is  (kg^u^  )"  =  (KU)"  (l^en^^"  *    ^'^'^ 
the  assumption  of  equilibrium  reads 

L-1/,.  ..\-l  ^.  ../idui 


(KU)-^(kg^u)-^  «  u/\^\ 


or 


The  quantity  u  in  (2)  is  the  result  of  an  averaging  operation; 
consider  this  average  to  be  an  ensemble  average;  it  is  reasonable, 
and  consistent  with  experience  in  classical  statistical  mechanics, 
to  assxame  that  if  (2)  holds  on  the  average,  it  holds  for  most 
systems  in  the  ensemble,  i.e.,  for  most  flows  there  exists  a  range 
of  k's  typified  by  k^  such  that 


(KU)-^I^I  «  kg^u^  ,   u  -j/u^dx  . 

This  further  assumption  will  be  at  least  partially  justified  below. 

The  purpose  of  the  present  lecture  is  to  derive  a  number  of 
conclusions  from  the  assumption  (2),  used  in  conjunction  with  the 
equations  of  motion.   In  particular,  the  form  of  the  inertial 
range  spectrum,  i.e.  the  equilibrium  range  spectrum  in  the  limit 
of  very  small  viscosity,  will  be  derived,  and  a  physical 
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explanation  of  the  reasons  for  its  existence  will  be  presented. 
The  inertial  range  spectra  will  be  studied  in  three  cases:  The 
one-dimensional  Burgers  equation,  and  the  Navier-Stokes  equations 
in  two  as  well  as  three  space  dimensions. 

Equilibrluin  and  inertial  ranges  of  the  Burgers  equation. 
Consider  an  ensemble  of  solutions  of  the  Burgers  equation 

(3)  u.  +uu^-  vu.^^  =  0  ,    -00  •<  X  <  +00  ,    t  >  0  , 

w     X     XX  — .    —  — 

with  u(x,  O)  in  C^.   v  is  a  r.mall  viscosity.   The  behavior  of  the 
Fourier  transform  u(k)  of  u{x),  in  the  limit  as  v  -♦  0  and  k  -*  oo 
(in  that  order)  can  be  readily  determined.  As  v  tends  to  zero, 
the  solutions  of  (3 )  tend  to  the  solutions  of  the  inviscid 
equation 

(4 )  u .  +  uu^  =  0 

in  L^;  u(k)  thus  tends  to  the  transform  of  the  solution  of  (4) 
pointwise,  and  uniformly  on  every  bounded  segment  of  the  k-£ixis. 
Hence  as  v  -►  0  the  behavior  of  u  is  determined  by  the  behavior  of 
the  solution  of  (4). 

The  solutions  of  (4)  exhibit  collections  of  shocks  separated 
by  intervals  of  smooth  u;  their  transforms  are  thus  of  the  form 

,  ia.(t)   , 

(5)  o(k-^)  +^  C^.(t)e  ^       k"-^  ; 

the  C^,  a^  are  functions  of  t  and  depend  on  the  initial  data; 
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they  are  in  this  sense  random.  The  form  (5)  is  universal,  i.e.   .^, 
independent  of  the  particular  initial  data.  For  large  k  the 
energy  spectrum  E(k)  is  thus  of  the  form        .  ^ 

(6)  ,^  E(k)  =  0(k-^)  , 

which  is  a  well  known  result. 

Let  us  now  rederive  (5)  and  (6)  using  the  assumption  (2). 
Take  the  Fourier  transform  of  (3).   We  already  know  that  if  we   are 
interested  in  the  inertial  range,  i.e.  in  the  equilibrium  range  in 
the  limit  of  decreasing  viscosity,  we  may  set  v=  0  in  the  resulting 
equation,  obtaining 

(7)  u^+ ikjru(k')u(k-k')dk'  =  0  , 

where  u(k)  is  the  complex  conjugate  of  u(-k). 

Let  k  be  in  the  inertial  range,  i.e.  of  order  k^  ,  and  make 
the  change  of  variable 

(8)  a*  =  O/U  ,   k*  =  k/K  , 

U,  K  defined  in  the  introduction.   Substitution  into  (7)  yields 

(9)  (KU)"^u*  =  ik*  J^u*(k')u*(k-k')dk'  . 

}'  ' 

p 
The  right-hand  side  is  of  order  k^  u  ,  the  left-hand  side  of  order 

(KU)"-'-  -g^  ;  by  (2),  for  most  flows 

(10)  ik  rG*(k«)u*(k-k')dk'  =  0  . 


Ill 
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The  energy-containing  eddies  have  been  relegated  to  the  neighbor- 
hood of  the  origin,  and  the  integral  must  in  general  be  interpreted 
in  a  principal  part  sense.   In  the  unstarred  variables  (10)  reads 


(11)  lim  ik  /  u(k-k«  )u(k')clK:'  ■=  0  . 


k~*oo 

Equation  (lO)  is  readily  solved:   u*  is  the  traneform  of  u*,  which 
satisfies 

(u*^)^  »  0 

u*^  m   c*^  ,   C*  =  constant, 
or 

u*  -  ±C* 

where  different  signs  may  be  assumed  on  distinct  parts  of  the 
x-ajcis.   Only  a  solution  with  a  single  change  of  sign  can  satisfy 
the  entropy  condition;  its  Fourier  transform  is  of  the  form 

(12)  G*  -  Ce^^^k-^  , 

where  C,  a  may  depend  on  t.  We  now  observe  that  any  super- 
position of  solutions  of  the  form  (12)  satisfies  (11);  Indeed  let 

ia  k  _,      ia  k 
u  'v,  C,e     k   +  C^e      as  k  -♦  oo  , 


Then 
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P  ia,k  ia^k      n   ,        -, 

j  G(k')a(k-k')dk«  =  (Cj_e   ^  +  c^e      )  J  p- ra^ '*^' 

ia^k    r    i(a,-a   )k'     ,       3^ 

Perforin  the   change   of  variables  k'    =  Kk* ' ,    k  =  Kk*;    we  obtain 

J  u(k')u(k.k')dk'   =^[(c^e     1     +0^6  )jl^pr^dk*' 

la  k    r     i(a  -a    )Kk*'       ,  ^  -1 

as  k-^oo,  K-»oo,  the  last  integral  is  negligible  compared  to  the 
first,  and  (11)  is  satisfied.   Thus  the  solution  of  (11)  is 

1  la.(t)   , 

o(k"  )  +  5~  C  (t)e  J    k~   as  k  —  00  , 

in  agreement  with  (5).   By  comparison  of  two  arguments  which  lead 
to  (5)  we  see  that  the  equilibrium  hypothesis  holds  for  all  flows; 
we  have  furthermore  an  interpretation  of  the  equilibriiim 
hypothesis:   As  the  flow  evolves  and  the  values  of  u  on  each  side 
of  a  shock  vary,  the  shock  adjusts  to  the  change  instantaneously. 
It  is  readily  seen  that  (5)  represents  an  equilibrium  in  wave 
number  space.  '■^ 

Two  dimensional  Navier-Stokes  equations.   It  has  been  shown 

2 
in  some  generality  (see  Ebin  and  Marsden  and  the  references 

therein)  that  as  the  viscosity  v  tends  to  zero  in  the  absence  of 
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boundaries  the  solutions  of  the  Navier-Stokes  equations  with 
smooth  initial  data  tend  to  the  corresponding  solutions  of  the 
Euler  equations  strongly,  as  well  as  in  L^^.  The  respective 
Fourier  transforms  then  tend  to  each  other  uniformly  on  every 
bounded  region  in  k-space  and  the  limit  v  ->  0  may  be  studied  by 
setting  V  =  0.   One  then  performs  the  scaling  analogous  to  (1) 

(13)  k*  «  k/K  ,   u*  =  u/U  , 

appeals  to  the  equilibrium  hypothesis  (2),  and  obtains  the  steady 
inviscid  equation 


where 


Qg^  «J'a*(k-k')u*(k')dk'  , 


k  k 

V/  •=  ^ay   "  -^  '      ^^ay  ^^^   Kronecker  delta), 

and  the  summation  convention  is  in  use.  u*(k)  is  the  complex 
conjugate  of  u*(-k),  and  the  pressure  has  been  eliminated  through 
the  use  of  the  equation  of  continuity 


(15)  k„a;  =  0  , 

giving  rise,  in  the  well  known  manner,  to  the  projection  P  . 

ay 

The  integrals  Qg  must  of  course  be  interpreted  in  a  principal-part 
sense.  As  before,  (l4)  can  be  solved  by  application  of  the  Fourier 
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transform:   it  is  the  transform  of  the  steady  Euler  equations, 
which  admit  an  infinite  number  of  solutions,  in  particular  circular 
vortices  whose  stream  function  has  a  Fourier  transform  of  the  form 

I,   II    II 

(16)  a*  =  -ik^/k^^  ,    u*  =  ik^/k^^  ,    k  =  /kf+k|  ,  P  constant; 

given  p,  the  energy  spectrum  is 

E(k)  =  0(k''^)  ,    o  =  4p-  3  . 

These  vortices  play  the  role  played  by  the  shocks  in  the  one 
dimensional  problem.   A  lower  bound  for  a  can  be  readily  found: 
The  H,  norm  of  u  is  constant  (by  conservation  of  vorticity),  thus 

k^E(k)dk 


/ 


is  uniformly  bounded,  and  a  >  3.  The  following  conjecture  appears 
natural:  in  any  particular  flow  vortices  will  appear  with  spectra 
of  the  form 

E(k)  =  0(k'^"^)  ,   e  -  e(t)  >  0  .  ,  -. 

When  E(k)  decays  as  fast  as  k''  or  faster,  finite-difference 
computations  have  a  chance  of  yielding  reasonable  results;  such 
computations  have  been  carried  out,  e.g.  by  Zabuski  and  Deem  , 
and  lend  support  to  the  conjecture.   In  an  ensemble  of  flows,  the 
behavior  of  E(k)  for  large  k  will  be  determined  by  the  vortices 
whose  transforms  decay  most  slowly,  and  should  then  be 
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indistinguishable  from  0(k''^),  as  postulated  by  Leith  and 
Kraichnan^ . 

The  three  dimensional  Navier-Stokes  equations.   In  three 
space  dimensions,  the  existence  of  solutions  of  the  Navier-Stokes 
and  Euler  equations  for  all  time,  and  a  fortiori,  the  convergence 
of  the  former  to  the  latter,  have  not  been  established.  We  have 
to  assume  that  the  solutions  exist  and  that  such  convergence  takes 
place.   These  assumptions,  together  with  the  equilibrium  assump- 
tion (2),  lead  to  the  three-dimensional  analogue  of  (l4). 

It  can  be  shown  that  if  grad  u*  ?^  0,  equations  (l4)  in 
three  dimensions  admit  only  two  dimensional  solutions.   It  should 
be  remembered  that  u*  characterizes  the  behavior  of  u  for  large  k. 
The  fact  that  u*  is  two-dimensional  does  not  imply  that  u  (and  u) 
are  two-dimensional.   The  u*  are  vortices;  their  "two- 
dimensionality"  means  that  their  radius  of  curvature  is  of  order 
^en  ^^^   that  the  flow  parallel  to  their  axis  is  relatively  smooth. 
Vortex-stretching  is  thus  allowed. 

We  consider  now  solutions  of  the  form  (l6).   The  non- 
increasing  character  of  the  energy  implies  that  a  >  1.   If  it 
were  known  that  the  integral 


J      k'^E(k)( 


remains  tiniformly  bounded  for  a  ■<  a  while  it  may  become  unbounded 
for  a  >  a^,  then  we  could  conclude  that  a  '=  a^+  1.     Such  results 
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are  unforttinately  unavailable.  A  more  precise  determination  of  a 
requires  a  heuristic  argument,  in  the  absence  of  a  rigorous  one. 
Consider  the  question,  which  vortices  can  arise  from  smooth  initial 
data.   Clearly  only  vortices  which  are  in  some  sense  stable  can 
develop.  We  are  in  presence  of  an  unusual  stability  problem, 
u*  characterizes  the  behavior  of  u  for  large  k  only.  We  have  thus 
to  look  for  vortices  u*  with  the  following  property:  When  they 
are  imbedded  in  a  flow  u,  and  when  that  flow  is  perturbed,  the 
behavior  of  the  transforms  u  for  large  k,  i.e.  the  form  of  u*,  is 
unaffected.  The  shock  solutions  of  (10)  have  such  a  property:   if 
a  flow  with  shocks  in  one  dimension  is  perturbed,  it  remains  a 
flow  with  shocks,  with  transform  of  the  form  (5).   Perturbations 
containing  a  bo\inded  set  of  Fourier  components  are  of  no 
significsince  when  applied  to  equation  (14),  since  a  scaling  will 
relegate  them  to  a  neighborhood  of  k*  =  0  and  will  leave  (l4) 
unaffected.   To  formulate  a  reasonable  problem,  appeal  must  be  made 
to  physical  intuition.  A  flow  encountered  in  nature  will  have  a 
definite  scale,  based  either  on  the  scale  of  the  driving  mechanism 
or  in  the  size  of  the  container.  u(k)  will  then  have  for  support 
not  the  whole  k-space  but  only  a  lattice  of  wave  numbers  which  can 
be  excited.   Let  h  be  a  length  typical  of  the  lattice.  The  inte- 
grals in  (l4)  become  sums  over  the  lattice  points;  let  us  write 
the  resulting  equation  in  the  symbolic  form 


(17)  Fj^[u*]  -  0 
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Formally,  as  h -*  0,  (17)  converges  to  (1^).   In  analogy  with  (l6), 
we  seek  solutions  of  (17)  of  the  form 

^,1  =  -ik2/H(k,h)  ,   ^h,2  "  iiti/H(k,h)  , 

k  =  (k^^^kg)  ,   k  =  JkJ+  kg  , 

where  ^  is  defined  only  at  the  nodes  of  the  lattice,  and  H  is  a 
function  of  k  and  h.  To  each  h  will  correspond  a  spectrum  E(k, h); 
if  H  is  a  homogeneous  function  of  k  the  spectrum  will  be  of  the 
form 

E(k)  «  0(k'''^)  , 


where  a.    depends  on  h.  The  question  is:   as  h  -»•  0,  does  a.  tend 
to  a  definite  value  of  a,    independent  of  the  lattice?  The  smswer 
appears  to  be  affirmative,  with  <^u~*''^>    ^  "   1*7  ±0.1,  a  result 
compatible  with  the  Kolmogorov  law  (see  Chorin  ). 

Conclusion*  We  have  made  plausible  a  picture  of  turbulent 
flow  as  a  disorderly  process  in  which  small  scales  are  dominated 
by  a  tangle  of  vortices,  whose  presence  explains  the  existence  of 
an  inertial  range  and  the  form  of  the  corresponding  spectrum. 
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15-1 
Lecture  I5 

Dynamo  Theory 
Stephen  Childress 

Introduction 

In  1919  Larmor  put  forward  the  idea  that  magnetic  fields 
found  in  nature  may  in  some  cases  be  maintained  against  dissipa- 
tive  losses  by  magnetohydrodynamic  induction  within  an  icolated 
homogeneous  fluid  region,  in  much  the  same  way  that  a  magnetic 
field  is  generated  in  a  laboratory  dynamo.   In  particular,  it  is 
believed  that  the  earth's  magnetic  field  is  generated  and   main- 
tained in  this  way. 

The  dynamo  models  of  geomagnetism  may  be  distinguished 
mathematically  from  other  problems  of  magnetohydrodynamics  by  the 
fact  that  the  equations  and  boundary  conditions  permit  "trivial" 
solutions  with  zero  magnetic  field.   The  interesting  problems 
therefore  have  an  eigenvalue  or  bifurcation  character.  Also,  in 
the  absence  of  a  detailed  description  of  the  velocity  field  within 
the  earth's  fluid  core,  the  theory  is  in  some  ways  concerned  with 
the  inverse  problem  of  determining  the  structure  of  the  dynamo 
from  surface  measurements  of  the  resulting  field. 

Magnetic  surveys  over  the  last  I50  years  indicate  that  the 
surface  field  is  almost  completely  of  internal  origin,  and  that  it 
may  be  regarded  as  consisting  of  two  parts,  namely  a  main  dipole 
field  and  an  irregular  non-dipole  field.  From  studies  of 
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magnetization  in  sea  floor  sediments  and  volcanic  rocks  it  is  knovm 
that  the  main  field  has  reversed  in  polarity  about  200  times  over 
the  last  80  million  years.   Periods  of  tmiform  polarity  may  be 
divided  into  "events"  lasting  10  -10-^  years,  and  "epochs"  lasting 
10  years  or  so;  there  appears  to  be  no  preferred  polarity.  The 
time  scale  of  events  is  roughly  equal  to  that  of  the  free  decay 
of  the  magnetic  field  that  would  occur  if  the  dynamo  process  was 

not  occurring.  The  non-dipole  field  undergoes  "secular  varia- 

p  Li 
tions"  with  time  scales  from  10-10  years.  P.   portion  of  this 

can  be  interpreted  as  a  westward  drift  of  an  otherwise  "xigid" 

field,  at  a  rate  of  about  .2°  long. /year.  The  non-dipole  field 

has  in  some  studies  been  represented  by  8  or  10  regional  sources 

within  the  core,  convected  relative  to  the  mauitle  by  a  motion 

consistent  with  the  observed  westward  drift  velocity.   It  has  also 

been  suggested  that  the  secular  variation  includes  hydrodynamic 

waves  moving  relative  to  the  maantle. 

The  kinematic  problem 

Most  work  in  dynamo  theory  has  been  concerned  with  the 
kinematic  version  of  the  problem,  wherein  the  velocity  field 
u(x, t)  is  prescribed  throughout  the  core  region  D  and  the 
dynamical  equations  are  discarded.  The  rationale  here  is  that, 
if  u  generates  a  non-decaying  magnetic  field  h(x, t),  then  u,  h  is 
consistent  with  ajiy  dynamical  model  for  some  distribution  of 
externally  imposed  forces;  of  course  it  could  well  be  that  the 
required  forces  are  not  at  all  realistic,  so  the  usefulness  of  the 
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method  depends  on  a  shrewd  a  priori  choice  of  the  class  of  u  to 
be  considered. 

The  problem  is  defined  as  follows:  For  the  given  solenoidal 
u(x, t),  determine  real  functions  h(x, t)  and  E(x, t)  (electric  field) 
satisfying 

V  X  h  =  aE  +  ouu  X  h 

Vx  E  =  -u  1^  '  "^'^  ^   ° 

in  D,  {a,\x   «  electrical  conductivity  and  magnetic  permeability  of 
the  core  =  constants), 

V»h  =  Vx  h  =  V.E  =  0 

hh 


VxE  =  -p.  ^-T  >  r^{E|  and  r^|h|  bounded 


in  the  region  D'  exterior  to  D,  with  h  and  the  tangential  compo- 
nents of  E  continuous  on  hB,    such  that  the  total  magnetic  energy- 
exceeds  some  5  >  0  for  arbitrarily  large  t.   Presumably  u 
represents  a  source-free  motion  of  an  incompressible  fluid  with 
zero  flux  through  the  boundary. 

If  u,  h,  and  E  are  regarded  as  independent  of  time,  the 
equation  for  h  in  D  may  be  written 

V  h  +  oy.  Vx  (ux  h)  =  0  . 

In  this  case  the  kinematic  dynamo  problem  reduces  to  an  eigenvalue 
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problem,  with  a  typical  speed  of  the  dynamo  appearing  as  the  eigen  - 
parameter.  .    ■   '  .  -i   :   ,   ; . 

It  is  known  that  in  a  sphere  of  radius  L,  u  is  a  dynamo  only 
if  R  H  u^aixL  >  TT,  where  u  is  the  maximum  of  |u|  over  D  and  time. 
Apart  from  the  flexibility  obtained  from  parameters  implicit  in  u, 
this  inequality  limits  asymptotic  solutions  of  the  eigenvalue 
problem  to  the  case  of  large  magnetic  Reynolds  number  R.   Cowling 
(1934 )  and  others  have  shown  that  u  cannot  be  a  dynamo  if  the 
associated  h  is  axisymmetric.   Also  it  is  known  that  certain  simple 
choices  of  u,  for  example  a  purely  "toroidal"  motion  of  the  form 

u  =  7x  T(x)x 

csuinot  maintain  any  steady  field  (Bullard  and  Gellman  (195^  )•   See 
also  Roberts,  P.  H.  (1967)).   In  the  steady  case  Herzenberg  (195Q)» 
and  in  the  non-steady  case  Backus  (I958)  have  provided  examples  of 
kinematic  dynamos  in  a  spherical  core.   Backus'  model  consists  of 
three  time  intervals  of  rigorous  motion  separated  by  sufficiently 
long  periods  of  free  decay  of  h.   The  initial  field  is  decomposed 
into  a  principal  part  f  and  an  arbitrary,  bounded  residual  r.   The 
velocity  fields  are  then  chosen  so  that,  after  one  full  cycle,  the 
solution  of  the  dynamo  equations  has  a  similar  decomposition  f+r' 
where  r'  is  no  larger  than  r.   Herzenberg' s  model  consists  of  a 
finite  number  of  small,  rigid  rotating  spheres  imbedded  in  an  other- 
wise stationary  core.   For  small  ratio  of  sphere  to  core  radius  the 
eigenvalue  problem  can  be  solved.   In  either  case  the  analysis  is 
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based  upon  the  different  rates  of  decay  (with  respect  to  space  or 
time)  of  the  functions  used  to  represent  h. 

Another  approach  utilizes  a  complete  analysis  in  an  infinite 
core  of  induction  by  space-periodic  motions  (Childress  1970, 
Roberts,  G.  1970).  The  periodic  dynamos  obtained  from  such  an 
analysis  can  be  easily  imbedded  in  a  sphere,  smd  the  resulting 
steady-state  dynamo  problem  can  be  solved  by  perturbing  a 
comparison  eigenvalue  problem.  An  example  of  a  spherical  dynamo 
of  this  kind  is 

u  =  eVxco(  |x|  )v(x/e) 

where  cd  is  a  scalar  cut-off  and  v(x)  =  (sin  Xp  +  cos  x-,, 
sin  X,  +COS  Xp,  sin  x,  +cos  Xp).   The  parameter  e  is  the  ratio  of 
the  scale  of  the  periodic  field  to  the  radius  of  the  core.   If  e 
is  sufficiently  small,  such  a  motion  interacts  with  a  large  scale 
h  to  produce,  upon  averaging  over  the  small  scale,  a  mean  current 
roughly  parallel  to  h. 

Almost  symmetric  dynamos 

A  fonnal  asymptotic  procedure  for  R  >>  1  has  been  devised 

by  Braginskii  (1964a, b).   In  these  models  deviations  of  u,  h  from 

-1/2 
eixial  symmetry  are  of  order  R  '  ,  a  typical  choice  of  u  being  of 

the  form  (now  using  dimensionless  variables) 

u  =  W(z,r)i^  +R-^/^u'(z,r,(J>)  +  R"^x  [^(z,r)i^]  . 

Here  (z,  r,  <)))  denote  cylindrical  polar  coordinates  and  u'  has  zero 
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4-mean.   Braglnskii  shows  that  the  (j)-mean  of  the  dynamo  equations 
in  D  reduce,  to  leading  order  in  R,  to  the  following  system: 

^+r-\.V(rAj.v'^A^  =  FB  , 

||  +  rUg.v(r"S) -V^B  -  [v(r'^)xv(rAQ)}^  . 


He 


re  the  (|>-mean  of  h  is  equal  to  h^  =  B(z,  r)i  i  +R"^  x  [A(z,  r)i  i  ], 

i/  -f  A  -A 

u  =  Vx  [i^g(z,  r)ii],  and  r,  -4* — ,  and  — g—  are  each  i^-means  of  a 

function  quadratic  in  u' . 

Braginskii's  equations  are  remarkable  in  that,  if  the  term 
involving  r  is  deleted,  the  resulting  system  is  identical  to  the 
equations  for  an  axially-symmetric  dynamo  involving  the  "effective" 
quantities  ^  and  A  .   The  surviving  term  involving  the  toroidal 
velocity  field  W  acts  as  a  source  of  the  azimuthal  component  B  of 
h^.   (This  "W-effect",  which  can  be  interpreted  as  a  drawing  out 
of  the  lines  of  force  of  the  meridional  component  Vx  (A^ii ),  is 
likely  to  be  an  important  one  in  a  rotating  spherical  core,  since 
\J   will  then  contain  the  geostrophic  component  of  velocity. )  The 
new  "r-effect"  provides  a  corresponding  source  of  meridional 
field,  thus  completing  the  cycle.   Numerical  solution  of  the 
eigenvalue  problem  for  steady  dynamos  has  been  carried  out  for 
several  choices  of  W  and  ^  .   In  certain  cases  the  ratio  of  the 
non-axisymmetric  to  axisymmetric  part  of  the  exterior  meridional 
field  is  found  to  be  of  order  R"  '  ,  in  which  case  R"  '  can  be 
interpreted  as  the  order  of  magnitude  of  the  tilt  of  the  geo- 
magnetic dipole  relative  to  the  earth's  axis  of  rotation. 
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Lecture  l6         ,   •   >:,        '". 

Linear  and  Nonlinear  Magneto-Elastic  Wave  Motion 

.':     ■  I    ■'•■!"  J,  Bazer  "   :■,■•''■    ■'<''■       .       •        'i        '■''' 

The  following  is  an  abstract  of  a  talk  on   magneto-elastic 
wave  motion  based  on  the  papers  [l]-[3]  listed  below  in  the   '  -v 
reference  section. 

An  account  is  given  of  some  work  dealing  with  wave  motion, 
linear  and  nonlinear,  in  an  electrically  neutral,  infinitely- 
conducting,  elastic  medium  subject  to  magnetic  as  well  as  elastic 
stresses-   The  analysis  is  based  on  the  fact  that  the  governing 
equations  of  such  a  physical  system  share  with  the  equations  of 
conventional  elasticity,  gasdynamics  and  magne togas dynamics 
the  property  of  being  representable  as  a  first-order  symmetric- 
hyperbolic  system  of  conservation  laws.  This  property  implies,  in 
advance,  that  the  linearized  magne to- elastic  equations  possess  a 
geometrical  theory  analogous  to  geometrical  optics  and  that  the 
nonlinear  equations  possess  a  theory  of  shocks  and  simple  waves  ' 
similar  to  that  of  gasdynamics  and  magne togas dynamics. 

A  complete  solution  of  the  following  basic  problem  in  the 
geometrical  theory  is  presented:   Let  vi   denote  a  surface      ''• 
bounding  a  small  disturbance  in  an  otherwise  undisturbed  medium. 
Suppose  the  discontinuity  across  5/  of  the  meignetic,  velocity 
and  strain  fields  are  specified.  The  problem  is  to  determine  the 
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wave  fronts  which  evolve  frcan  sc^,^  —   there  will  be  six  of  them  in 
general  —  and  to  ascertain  the  orientations  of  the  field  quantities 
and  the  strengths  of  the  discontinuities  carried  on  these  fronts. 
Insofar  as  the  propagating  fronts  are  concerned,  the  presence  of 
the  magnetic  field  is  reflected  in  the  fact  that  their  propagation 
speeds  are  direction-dependent.  Specifically,  these  speeds  depend 
upon  the  angle  between  the  direction  of  propagation  and  that  of 
the  local  magnetic  field.  As  to  the  orientation  of  the  field 
quantities  on  the  fronts,  these  are  found  to  be  unambiguously 
determined  by  the  direction  of  propagation  and  the  local  magnetic 
field,  provided  these  are  not  parallel.  Moreover,  in  only  one  of 
the  three  waves  that  may  propagate  in  a  given  direction,  the  so- 
called  intermediate  wave,  are  the  velocity  and  displacement 
vectors  purely  transverse,  that  is,  parallel  to  the  wave  front. 
In  the  remaining  'slow'  and  'fast'  fronts  both  longitudinal  and 
transverse  components  are  present.  These  results  are  completely 
analogous  to  those  of  geometrical  hydromagnetics  and  in  sharp 
contrast  to  those  of  conventional  geometrical  elasticity  both  of 
which  incidentally  are  limiting  cases  of  the  present  theory. 

As  in  gas  dynamics  the  nonlinear  theory  is  restricted  to 
one- dimensional  wave  motion.  This  means  that  all  physical 
quemtities  are  required  to  depend  on  only  one  space  variable  etnd 
time.  No  restriction  is  placed  on  the  orientations  of  the 
magnetic,  velocity  and  strain  fields.   In  the  simplest  case 
considered,  the  stress-strain  relation  is  assumed  to  be  that  given 
as  Hooke's  law.  The  nonlinearity  is  then  due  to  the  interaction 
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of  the  magnetic  field  with  itself  and  with  the  velocity  and  strain 
fields.   Despite  the  complexity  of  the  equations,  it  is  possible 
to  obtain  almost  completely  explicit  simple  wave  solutions.  As 
in  raagnetogaedynamics  it  turns  out  there  are  slow,  fast  and 
Alfven-like  simple  waves  which  together  with  the  corresponding 
magneto- elastic  shocks  enable  one  to  solve  one-dimensional 
propagation  problems  with  sufficiently  simple  initial  and  boundary 
conditions.   One  interesting  result  of  the  nonlinear  theory  is 
that  magneto-elastic  simple  wave  motion  affords  a  mechanism  for 
generating  intense  maignetic  fields. 
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Lecture  17 

Viscosity  of  the  Earth's  Crust 
Gerald  Grube 

Introduction 

Geophysicists  have  observed  that  for  forces  of  long  duration 
the  earth  behaves  like  a  fluid  with  a  very  high  viscosity.   This 
phenomena  is  not  surprising  when  we  notice  that  ice,  pitch,  and 
glass  have  the  same  property.   Viscosities  for  many  substances 
commonly  thought  of  as  solids  may  be  found  in  various  handbooks. 

During  the  last  ice  age,  vast  ice  sheets  covered  large 
portions  of  the  earth  for  an  extended  period  of  time.   In  1928, 
Nansen  published  his  detailed  study  of  the  Scandanavian  ice  sheets 
and  the  resulting  deformation  of  the  earth's  surface.   Haskell 
derived  a  formal  solution  for  a  viscous  fluid  in  1935  and,  using 
Hansen's  data,  determined  a  kinematic  viscosity  for  the  earth's 
crust  of  5x  10   cm  /sec. 

Until  recently,  only  a  few  other  attempts  were  made  to 
calculate  a  viscosity  for  the  earth.   Haskell's  results  have  been 
used  by  many  geophysicists  in  their  theories  of  the  earth's 
interior.   Particularly,  the  value  of  viscosity  is  critical  to 
understanding  the  continental  drift.   In  the  following,  Haskell's 
work  will  be  summarized  and  a  discussion  of  how  one  might  try  to 
improve  upon  his  results  will  be  presented. 
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Formulation  ■': 

We  will  consider  a  semi-infinite,  incompressible  highly 
viscous  fluid  with  the  z  =  0  plane  as  the  equilibrium  surface. 
Cylindrical  coordinates  will  be  used  and  only  circularly 
symmetrical   distiiKbances  will  be  considered.   Because  the 
viscosity  v  is  large,  we  may  neglect  inertial  and  nonlinear  terms; 
the  equations  of  motion  are: 

fi  a 


r  er  (-  ^  ^r)  -  -^  ^r  ^  !^  \ 

r       oz 


=  -^  P  , 
dr   ' 


(1) 


r  or  ^  Cr     z 


^       ^2 

g2 


^ 


p  > 


i  ^  ( rV  )  +  V  =0 
r  cr  ^   r     z 


where  V  and  V  are  radial  and  horizontal  velocity,  respectively, 
P  =  p~  P  -  gz  with  p  the  pressure,  p  density,  and  g  the  gravita- 
tional force.   The  positive  z  axis  is  downward.  Assuming  small 
surface  deflections  4(r,  t),  vie   may  linearize  the  usual  boundary 
conditions  of  zero  stress  at  the  surface: 


S 


2v-^V^-P-g;;  =  0,    2  =  0, 


(2) 


5 


,   V  +  -^  V 
cr  z   cz  r 


0  , 


z  -  .,  , 


2 
^ 


z  ' 


z  =  0  . 


We  now  apply  the  Hankel  transform: 
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oo 

fJ^(^,z,  t)  =  J  rJ^(Ar)f(r,z,t)dr 


0 

With  n  =  1  for  V  and  n  =  0  for  the  other  variables,  one  may 

r 

transform  (l)  and  (2)  and  solve  the  resulting  problem  directly. 
In  particular  one  obtains: 

00 

V2(r,z,  t)  =y  K(A)(l+Az)e-''^-S*/^^^  J^(Ar)dA  , 
e(r,t)  =f  /  >.K(>)e-^^/^^^  J^Ur)dX 


0 


where  K{'K)   is  as  yet  undetermined.   Notice  that  K  is  known  if 
either  4(r,  O)  or  V^(r,  0,  O)  is  given. 


Geological  Data 

Nansen  has  determined  past  surface  elevations  and  rate  of 
uplift  at  two  sites,  Oslo  and  Angermanland,  about  50O  km. 
apart;  the  latter  having  been  near  the  center  of  the  ice  sheet. 
We  assume  a  reasonable  surface  form: 


(4)  e(r,0)  =  :3e-^  ^^1-bV)  . 


Substituting  (4)  into  (3),  one  may  solve  for  v: 


<5)  V  =  -.22  bV,(0',0,O)  * 

The  parameters  b  and  3  are  determined  so  that  (4)  fits  the  surface 
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elevations  given  by  Nansen,  then  the  rate  of  uplift  determines  v. 

Using  data  from  various  times  in  the  past,  Haskell  obtains 

21    2 
V  =  2.9x  10   cm  /sec  as  an  average  value. 


Acceleration  Terms 

One  may  justify  not  including  the  acceleration  terms  in  (l) 
by  means  of  a  singular  perturbation  (a  new  unstretched  time  is 
defined  by  multiplying  the  old  time  by  a  large  parameter).   Of 
course,  now  one  is  restricted  in  terms  of  initial  conditions; 
either  |(r, O)  or  V  (r, 0, O)  completely  determines  the  solution. 

However,  one  may  include  the  time  derivative  terms  in  (l) 
and  still  obtain  a  closed  form  for  the  solution.  We  proceed  as 
before  but  now  the  Laplace  as  well  as  the  Hankel  transform  must 
be  applied: 

GO  oo 

f^(A,z,s)  =  J  ^Jn^^""^/  e"^^f(r,  z,t)dtdr  . 
0  0 

We  use  Haskell's  solutions  (3)  to  provide  initial  values  for  V^ . 

Then,  the  solution  becomes: 

CO      ioo 
V^(0,0,t)=^j   AdA  f   e^*  V^^^(A,0,s)ds 

0      -ioo 

(6)      ^,^,i^,0,s)    =  (1)^/2  K(A)(i^-^) 


=  A+  (2a2+s)2-  4a5(a^+  s)^/2  . 
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The  function  A  is  common  in  viscous  problems  and  has  been  studied 
by  Lamb  (p.  627)  and  Miles.   By  bending  the  contour  of  integra- 
tion and  using  asymptotic  expressions  for  the  zeroes  of  A,  one 
can  simplify  (6)  into  a  form  more  suitable  for  numerical  evalua- 
tion on  a  computer. 

From  Nansen's  data,  additional  values  of  viscosity  may  be 

computed  by  studying  the  evolution  of  the  solution.   The  viscosi- 

21 
ties  obtained  in  this  manner  have  an  average  value  of  3.1x10 

2  / 
cm  /sec. 

Lake  Bonnevile 

The  difficulty  in  computing  a  viscosity  from  the  recovery 
of  the  earth  from  the  Scandinavian  ice  sheet  was  due  to  meager 
data.   The  flooding  of  the  Lake  Bonneville  basin  in  Utah,  however, 
was  an  event  for  which  we  have  some  idea  not  only  of  the  deflec- 
tions involved  but  also  the  history  of  the  surface  loading.   If 
one  adds  the  normal  stress  due  to  the  loading  in  (2a),  then  a 
solution  for  the  free  surface  similar  to  (6)  may  be  derived. 
Viscosities  may  now  be  computed  numerically  by  using  the  observed 

deflection  of  the  earth's  surface  near  the  lake  center  (see  Grube). 

20    2 
In  this  manner,  an  average  viscosity  of  3«6xlO   cm  /sec  was 

computed. 

Conclusions 

Although  Haskell's  approach  is  self-consistent,  it  is  a 
singular  perturbation  technique  and  all  the  initial  conditions 
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cannot  be  used.   When  the  acceleration  terms  are  included,  we 
still  may  derive  a  closed  form  for  the  solution,  which  may  easily 
be  evaluated  on  a  computer.   Thus,  if  one  possesses  enough  data, 
solutions  to  the  linearized  Navier-Stokes  equations  may  be  derived 
and  evaluated  numerically. 
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Lecture  18 

Nonexistence  of  Shear  Waves  in  the  Earth's  Core 
J.  J.  Stoker 

It  is  a  well-known  fact  of  seismology  that  the  core  of  the 
earth  does  not  transmit  shear  waves.   In  geophysics  this  is 
commonly  explained  by  postulating  that  the  core  of  the  earth 
behaves  like  a  liquid,  in  spite  of  the  fact  that  other  phenomena 
indicate  that  the  core  is  more  rigid  than  steel.   In  this  lecture  it 
will  be  shown  that  according  to  a  nonlinear  theory  of  elasticity, 
shear  waves  cannot  propagate  in  an  elastic  medium  if  the 
compression  is  sufficiently  large 

We  shall  consider  an  elastic  medium  with  plane  strain.   Let 
X,  y  be  the  Lagrange-coordinates  when  the  medium  is  unstrained, 
and  X,  Y  be  the  Euler-coordinates  when  the  medium  is  strained. 
The  local  deformation  is  described  by  the  Jacobian  matrix 


PMPi,)  = 


r  ^x  hx 

'Sx  '  15y 

hY  hY 

"Sx  '  3y 


When  the  distortion  is  a  rigid  body  displacement,  P*  P  (an  asterisk  de- 
notes the  transpose)  is  equal  to  1,    the  identity  matrix,  there  is  no 
strain.   According  to  linear  theory  of  elasticity,  when  the  distortion 
produces  small  extensions,  P*  P  is  approximately  _I  +  2E  +  ...  where 

E  is  the  stress  matrix,  given  by  E.  .  =  o  -v^  (^■5-^^  )  +  o  ■';^;;-  (x.-x.). 
—  ij   c;  ox .   11   d   GX^   J   J 

Thus,  in  the  nonlinear  theory  of  elasticity,  when  the  distortion 


136 


18-2 


produces  large  extensions,  P*  P  should  be  related  to  E  as 
E  =  f(P*'  P),  the  function  t{i)   should  have  the  properties  f(l)  -  0 
and  f'(l)  =  -i  so  the  nonlinear  theory  will  cover  the  linear  theory, 
¥e  shall  choose  f(0  =/C-l'   Thus,  we  assvime  the  following 
relationship  between  E  and  P        : 


E  =  /P*  P-I 


Note  that  P*  P  is  a  symmetric  matrix,  hence  there  exists  an 
orthogonal  matrix 

cos  9  sin  e 


C  = 


-sin  6  cos  9 


such  that  P*  P  =  (PC)  ,   Hence 


I+E  =  PC  . 

The  local  angle  of  rotation  9   is  given  by  tan  0  =  (^-- •^-)/(-^ +  -^) 
so  that  E  is  a  sjnnmetric  matrix. 

It  is  assumed  that  the  elastic  property  is  characterized  by  a 
strain  energy  density  function  W  such  that  the  potential  energy  of  the 
elastic  body  in  the  strained  state  is  //  Wdxdy.   The  Lagrange  stress 
matrix  Q  is  defined  by  Q^.  =  SVl/oP^  ..   Then,  using  the  principle 
of  virtual  work,  it  can  be  shown  that  the  strain  matrix  t_  is 
related  to  P  and  Q  by        . 

T  =  (det  P)"-*-  QP*  . 
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Newton's  second  law  of  motion  shows  that  in  the  absence  of  body 
forces  the  equations  of  motion  are 

S^X    ^^11  ^  ^^12 

p  ^  =  -ex-    3r-  ' 

h^Y  _   ^^21  ,  ^^22 

where  p  is  the  mass  density  in  the  unstrained  state. 

For  an  isotropic  medium  the  dependence  of  ¥  on  P  must  be  such 
that  \J   is  a  function  of  the  invariants  of  P  under  rotation.   Thus 

W  =  W(r,  s)  , 
where 

--iFi=[,i.|,^.(i-if)'f'\ 

Here  tr  and  det  refer  to  the  trace  and  the  determinant  of  a  matrix., 
Accordingly,  Q  and  t  are  given  by 


—   3r  —   ~  ds  — 


Q  =  ^.    C*+S  ^  P*-^- 


—   s  dr  •—  --•'   ds  — 

The  function  W(r,  s)  should  satisfy  the  condition  ^+4^  =  0  at 
r  =  2,  s  =  1  so  that  zero  stress  will  correspond  to  zero  strain. 
Now,  T  is  given  by 

(iM'i  -  i  1 H^^^  ^  i!s  4'i -i)  !'^  5'i  -  - 

\ or  os      / 
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when  W(r, s)  is  expanded  at  r  =  2,  s  =  1.   On  the  other  hand,  in 
the  linear  theory  ( P*  P  -^  1  hence  r  ->•  2,  s  -*  l),  t  is  equal  to 
A(tr  E)l+2nE,  where  A  =  Ev  (l+v  )"''■( l-2v  )""^  and  H  =  2  E(l+v)"   are 
the  two  Lame  constants  defined  with  respect  to  Young's  modulus  E 
and  Poisson's  ratio  v.   The  corresponding  function  W(r, s)  in  the 
limit  of  the  linear  theory  is  -^  (r-2)  +M-(r  -  2r-2s+  2),  which  is 

equal  to  ^  (tr  E)^  +|j.(E  E  ) .       '■   ''-   •'  "         ' ''  "' 

For  a  harmonic  medium,  as  defined  by  F.  John,  W(r, s)  has 
the  following  form   "  ,      ■  -i  --  ■• 

¥(r,  s)  =  2M.[F(r)  -  s]  . 

At  r  =  2,  the  function  F(r)  should  satisfy  the  conditions  F(2)  =  1, 

F'(2)  =  1,  F"(2)  =  (A+2|j.)/2u  in  view  of  the  fact  that 

F(r)  ->  1+  (r-2)  +  ^^i^  (r-2)'^+  .  . ,  in  the  limit  of  linear  theory. 

In  addition,  F(r)  should  satisfy  the  conditions  -j-  —  -r—  >  0  and 

d  P 

— ^  >  1,   The  first  condition  follows  from  the  fact  that  the  stress 

dr'^ 

should  increase  monotonically  with  the  strain  in  a  state  of  uniform 

stress.   The  second  condition  follows  from  the  fact  that  when  a 

long  rectangular  bar  is  stretched,  the  decrease  in  its  lateral 

thickness  is  slower  than  the  increase  in  its  axial  length.   In 

terms  of  F(r)  the  relations  for  Q  and  t  are   '       '  '  '   ''  '      ' 


Q  =  2M.F'(r)C*  -2nsP*""'" 
X  =  2n  I^IeI   (i+E)  -  2hI 


and  the  equations  of  motion  are 
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_0_  ^^X  ^  S^A  _  bB 
2\x   -^   ox  '  oy 

_D_  S^y  ^  SA  ^  SB 
2(1  ^^   ^   ^ 

where  A  =  F' (r)  cos  6   and  B  h  F'(r)  sin  6,    with  e(x,y)  the  angle 
introduced  earlier  with  reference  to  the  rotation  matrix  C^.  VJhen 
a  harmonic  medium  is  in  equilibrium,  A  and  B  clearly  satisfy  the 
Cauchy-Riemann  equations,  hence  they  are  a  pair  of  conjugate 
harmonic  functions  of  x  and  y.   It  folloi\rs  that  6   and 
■^  log  [F'(r)]  ,  equal  to  arcta,n  S  and  log  (A^+ B  )  -^^  respectively, 
are  another  pair  of  conjugate  harmonic  functions. 

There  is  in  general  one  value  r  of  r  at  which  F' (r)  =  0;  it 
corresponds  in  general  to  a  very  large  strain  —  of  the  order  -^-^ 
sa,y..   Since  A+  iB  is  an  analytic  function  of  x  +  iy  and  F'(r) 
vanishes  at  r  =  r  ,  the  zeros  of  F' (r)  in  the  body  must  be  isolated, 
unless  F'(r)  is  identically  zero;  it  then  follows  that  r  >  r 
everywhere  or  r  <  r  everywhere  inside  a  harmonic  median:  whj.ch  is 
in  equilibrium  without  body  forces. 

Now  consider  wave  motions  in  a  strained  elastic  medium  w:i  th 
X  equal  to  ax,  and  Y  equal  to  by,  correspondingly,  r  equal  to 
r  s  a+b,  and  0   equal  to  0.   Let 

X  =  ax  +  u(x,  y,  t)+  ...  , 

Y  =  by  +  v(x,  y,  t)  +  . ..  . 

The  equations  of  motion  linearized  v/ith  respect  to  the  strained 
equilibrium  state  are 
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S^v 


0   S  U  _  p„  /On  /  O  U   CI  V  \         _1_  p,  /  0>,  f  .^U   d 
r 

Plane  wave  solutions  u,  v  ~  exp  i(^x+r)y-ct)  exist  only  if  the 

wave  speed  c  is  equal  to  Crp  =  [  (2iJ./o  )F' (r°)/r°]  '      or  equal  to 

c^  =  [  (2(jl/o)F"  (r°)]  ^  ,   In  the  first  case,  the  displacement 

vector  (ujv)  is  perpendicular  to  the  wa.ve  normal  (^,r]);  it 

represents  a  transverse  (shear)  wave.   In  the  second  case^  the 

displacement  is  parallel  to  the  wave  normal;  it  represents  a 

longitudinal  (dilatation)  wave.   We  remark  that  <^l  ^  c^  because 

-3 -3—  >  Oo   If  the  initial  equilibrium  position  is  without 

dr  r  dr  i         j- 

strain  (a  =  1,  b  =  l),  the  two  v;ave  speeds  are  real  c     =  /\^/?, 
c,  =  y(A+2M.)/p;  these  are  the  classical  results  from  the  linear 
theory  of  elasticity.   Since  F"(r°)  is  always  positive,  longitudi 
nal  waves  can  always  propagate  v^hatever  the  initial  strain  is. 
However,  if  the  initial  state  of  equilibrium  has  a  sufficiently 
large  compression  strain  so  that  the  medium  is  in  a  subcritical 
state  r°  <  r  hence  F'  (r°)/r°  is  negative,  then  Crp  is  imaginary 
and  no  shear  waves  can  propagate  at  all. 
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Lecture  19 

Gravity  Waves  in  a  Thin  Sheet  of  Viscous  Fluid 
(work  of  C.  C.  Mei) 

E.  Isaacson 

Let  us  consider  a  two  dimensional  incompressible  viscous 
flow  down  a  plane  inclined  at  an  angle  0.      The  governing  equations 
are 

^u+^  V  =  0  , 


|^u+u|^u  +  v|pu+p-l|^p  =  g  Sin  e  +  p-V^ 


''     .^,u. 


s  5  ....  a     .   -1  a  „      „  .  -1  /a^    .  s^ 


^  v  +  u^  v+v^  v+p-  g^  P  =  -g  cos  e+p-^^l^  +  ^1  V  . 

Here  u  is  the  velocity  in  the  x-direction  along  the  inclined  plane, 
V  is  the  velocity  in  the  y-direction  normal  to  the  plane,  p  is  the 
pressure.   The  constants  p,  g,  and  m-  are  the  density,  the  gravita- 
tional acceleration  and  the  coefficient  of  viscosity.   The  boundary 
conditions  to  be  satisfied  are 

u«=0,   v  =  0,on  the  inclined  plane   y  =  0  ; 

while  on  the  free  surface,  y  =  ri(x,  t),  we  have  the  kinematic 
condition, 

^n+u^Ti  =  v,on   y  =  Ti(x,  t)  , 
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and  the  condition  of  no  stress,  which  for  the  two  components 
yields,  -■    . 


(p-2u  ^  u)  ^  Ti+u(^  u  +  J^  v)  =  0  ,   on  y  =  Ti(x,t)  , 


i-  u)  ^  T1+|I(3-  u  +  ^^ 
dx  '  dx  '   ^Sy    6 


p-  2u  ^  v  +  u(^  u  +  ^  ^^  ^  "1  =  0   '   on  y  =  Ti(x,t)  . 
The  above  set  of  equations  admits  the  following  solution: 
Ti  =  h  , 

u  =  U(y)  =  (x'^g  sin  6   (hy  -|  y^)  , 
V  =  0  , 
P  =  P(y)  -   pg  CCS  e    (h-y)  , 

which  represents  a  steady,  stratified  flow  of  constant  depth. 

To  consider  a  nonlinear  time- dependent  flow  superposed  on 
this  steady  flow,  we  let 

u  =  U(y)  +u'  , 

P  e  P(y)  +  p'  . 

Now  consider  a  flow  with  time  scale  T,  length  scales  L  and  H  in 
the  X-  and  y-directions  (with  U.  =  h),   velocity  scale  V  and  such 
that  e  =   H/L  «  1  (shallow  fluid),  the  Reynolds  number  pu"  LV  =  e, 
the  Froude  number  V(gL)''-^/^  =  t^^^   (whence  V  =  (gH)"'"''^^),  and  the 
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-]_  -1     2 
Strouhal  number  LT  v   =  e  •   It  can  be  shown  that  this  implies 

2     1/2 
that  the  steady  flow  is  slow,  i.e.  U(y)  =  0(e  (gH)  '^  ),  and  that 

^  u'  and  ^  v'  are  O(e^),  u  ^  u'  +  v»  ^  u  and  u  ^  v+  v'  ^  v' 

are  0(e),  ^  p',  •$-  p',  ixV  u'  and  m-V  v'  are  0(e°).  Accordir^ly, 
it  can  be  found  by  a  formal  expansion  procedure  first  in  powers  of 
the  variable  y  and  then  in  powers  of  the  parameter  e,    that 

For  a  progressing  wave  of  permanent  form,  in  which  rj  is  a  function 
of  4  =  x-Ct,  T^  satisfies  the  following  equation 

^  (-Ce^Ti+^  sin  e   Tp-^  cos  e   rp  ^)   =  0   . 

Suppose   that  drj/d^  vanishes  at  t]  =  e,    then 

cos  0  Tj-^  ^  =  sin  e    (rp  -  e^)  -3Ce^(ri-e)    . 

Since  drj/d^  is  negative  when  r]  is  between  e  and 
b  =  ^  [-l+(l2C/sin  0-3)    ],    t]   decreases  montonically  between 
these  two  values.   In  the  critical  case,  b  =  0  (which  occurs  for 
C  =  -c-  sin  0),  the  free  surface  is  normal  to  the  bottom  (i.e., 
dri/d^  becomes  infinite  at  r)  =  O).   This  solution  represents  a  bore 
invading  a  dry  bed.   The  equations  described  above,  govern  the 
"large"  perturbations  about  the  uniform  flow.   Such  a  theory  may  be 
of  interest  in  connection  with  the  flow  of  glaciers.   Other, 
"smaller"  perturbation  equations  and  progressing  wave  solutions  are 
also  studied  in  C.  C.  Mei,  (1966),  Journal  of  Mathematics  and 
Physics,  45,  pp.  266-288 « 
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Lecture  20 

Temperature  Profile  of  the  Solar  Wind  j 

Tyan  Yeh 

The  solar  wind  is  the  most  important  discovery  in  space 
physics  in  the  last  decade.   It  refers  to  charged  particles 
moving  quickly  away  from  the  sun  in  the  interplanetary  space. 
This  solar  corpuscular  radiation  is  emitted  from  the  surface  of 
the  sun  due  to  its  high  temperature  (about  2x10  K  )  in  all 
directions  and  at  all  times.   The  operating  mechanism  of  heating 
in  the  outer  corona  is  so  effective  that  the  high  temperature  of 
the  corona  extends  well  into  interplanetary  space.   The  particles 
in  the  solar  wind  are  accelerated  in  their  flight  by  the  pressure 
gradient  just  like  a  gas  discharged  into  a  low  pressure  region 
through  a  de  Laval  nozzle.   At  the  orbit  of  the  earth  the  solar 
wind  has  a  flow  velocity  of  about  300  km/sec  with  a  particle 
density  of  about  5  ions/cm   and  a  temperature  of  about  10  K. 
The  interaction  of  the  supersonic  solar  wind  with  the  geomagnetic 
field  produces  a  magnetosphere.   Thus  the  dipole  geomagnetic 
field  is  distorted  considerably  into  a  configuration  with  a 
magnetic  tail,  and  the  geomagnetic  field  is  confined  inside  the 
magnetosphere. 

The  characteristic  supersonic  expansion  of  the  solar  wind 
was  predicted  successfully  by  Parker's  theory  of  coronal  expan- 
sion.  Space-probe  measurements  confirmed  that  the  solar  corpus- 
cular radiation  forms  a  wind  rather  than  a  breeze  as  proposed  in 
the  evaporation  theory.   In  the  hydrodynamic  theory  of  coronal 
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expansion,  the  velocity  profile  of  the  solar  wind  is  represented 
by  a  critical  solution  of  the  momentum  equation,  assuming  that 
the  temperature  is  known.   The  temperature  profile  itself  is 
determined  by  the  energy  equation.   Beyond  several  solar  radii, 
thermal  conduction  is  effectively  the  sole  means  of  transport  of 
heat.   Thus  the  governing  equations  for  a  steady  state,  spheri- 
cally symmetric  solar  wind  are 

(1)  ^  (nur^)  =  0  , 

(2)  (mp+mg)nu  -ff  +  ^  (2kTn)+  (m^+mg  )GM  ^  =  0  , 

-,   ,  /  p  m  +m_   p        \   1   H    Q 
^  ^  '  ^  (-^^  nu^+3kTn)u)+-^  ^  (r^(2kTn)u)  +  (m  +m^  )GM  ^n 

(5)  :    , 

_  1   d   /„2   5/2  dT^ 

-  -^  d?  ^"^  ^^     d?^  • 
r 

Here  r  is  the  distance  from  the  center  of  the  sun,  n  is  the  number 
density  of  protons  or  electrons,  u  is  the  flow  velocity,  T  is  the 
temperature.   The  constants  m  ,  m   and  M  are  the  masses  of  a 
proton,  an  electron  and  the  sun,  k  is  the  Boltzmann  constant,  and 

G  is  the  gravitational  constant.   The  thermal  conductivity  of  a 

S/2 
fully  ionized  gas  depends  on  the  temperature  k   =   KT  '  ,  K  being  a 

constant.   The  pressure  is  equal  to  2kTn,  the  internal  energy  is 

equal  to  3kTn. 

The  continuity  equation  (l)  can  be  integrated  to  give 

(4)  nur  =  F  . 
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The  constant  4TrF  is  the  total  particle  flux.   Then  the  momentum  "*■ 
equation  (2)  can  be  written       '■  '  '■  '■"■ 

.    ...  •  '      j:-i  '  r-}'' 

'  r  du        2  -  (r/T)(dT/dr)  -  (mp+mg)GM/2kTr 

^  '^'^  (m  +mg)uV2kT-l         * 

-z 

The  left  side  of  equation  (3)  is  equal  to  2k(^  nu  dT/dr  -  Tu  dn/dr) 
by  virtue  of  equations  (l)  and  (2)o   Thus  the  heat-flow  equation 
(3 )  can  be  written 

(6)  JL^  fr^T^/^  d^)  _  m  _  T  ,    r  dus 

lo;  2kF  dr  ^        dr'^    2  dr    r  ^    u  dr'' 

With  the  temperature  profile  assumed  known,  equation  (5)  was 
used  by  Parker  to  demonstrate  the  supersonic  expansion  of  the  solar 
wind.   The  coupled  equations  (5)  and  (6)  for  u  and  T  have  not  yet 
been  completely  solved  analytically.   So  far  only  special  solutions 
for  some  values  of  the  parameters  have  been  obtained  by  numerical 
integrations.   In  carrying  out  the  numerical  methods,  difficulties 
in  computation  do  not  arise  from  the  coupling  between  the  energy 
equation  and  the  momentum  equation;  they  arise  rather  from  the 
singular  boundary  condition  at  infinity  for  the  energy  equation. 
These  numerical  solutions  are  isolated,  because  it  is  not  clear  at 
all  how  they  are  imbedded  in  all  possible  solutions.   However  by 
assuming  that  the  velocity  profile  is  known,  equation  (6)  can  be 
used  to  study  the  temperature  variation  in  the  solar  wind. 

When  u  is  regarded  as  a  known  function  of  r,    equation  (6) 
is  a  second  order  differential  equation  for  T,  and  ru~   du/dr  as 


147 


c:0-4 


given  by  eqaation  (5)  contains  no  second  uerivative  of  T.   Thus, 
the  mathematical  structure  of  equation  (6)  is  effectively  inde- 
pendent of  the  variation  of  the  term  ru"  du/dr,  particularly  so 
when  du/dr  is  everywhere  positive  for  the  solar  wind.  We  shall 
exploit  this  mathematical  remark  by  assuming  that  ru"   du/dr 
is  equal  to  a  constant.  The  value  of  this  positive  constant  does 
not  affect  the  solutions  qualitatively.  Hence  we  shall  set  it  tc 
zero.  Formally  such  an  approximation  is  justifiable  when  the 
right  side  of  equation  (5)  is  small,  namely,  if  the  solar  wind 
has  a  small  logarithmic  expansion  rate.  At  large  distamces,  where 
the  solar  wind  attains  a  large  velocity  and  a  small  temperature, 
this  approximation  is  certainly  valid. 

W3  proceed  to  solve  the  approximate  heat-flow  equation 

/7^  K   d  /„2t5/2  dT^  3  dT  o  T  _ 

Two  boundary  conditions  are  required  to  fix  the  solution.   One  of 

them  arises  naturally,  viz.,  T  vanishes  at  infinity.  The  other 

one  can  be  chosen  as  prescribing  the  temperature  at  some  distance, 

say  T  =  T^  at  r  =  r  . 
"      o        o 

A  direct  verification  will  show  that  (2K/55kF)rT^^  =  1   is 
an  exact  solution  of  equation  (7)«   This  isolated  solution  is  use- 
less, because  it  cannot  satisfy  the  arbitrarily  prescribed  tempera.- 
ture  at  some  distance,  although  it  satisfies  the  sing-.rlar  condition 
at  infinity.   However,  upon  introduction  of  the  following  new 
variables 
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[Q)  X  =  j^^   rT 

(9)  y^.^r^gll, 

the  second  order  differential  equation  (7)  is  split  into  two 
simultameous  first  order  differential  equations: 

(10)  r^  =  x-y, 

nn  dy  _  14 y^  +  6y  -  20x 

(11)  r^ ^5^3j . 

The  advantageous  feature  of  the  simultaneous  equations  is  that 
elimination  of  r  results  in  a  single  first  order  equation 

f^o^                                          ^  -   1%^+  6y-  20x 
^^^^  c[x 35x(x-yJ —  • 

And  according  to  definitions  (8)  and  (9)  we  have 

f-i-zs  r  dT  ^   2  y 

Thus,  after  solving  equation  (12)  to  obtain  y  in  terms  of  x,  we 
can  solve  equation  (lO)  to  obtain  x  in  terms  of  r,  then  T  can  be 
found  in  terms  of  r  by  solving  equation  (13 )•   In  this  way.  all 
possible  solutions  of  equation  (7)  are  found;  and  the  particular 
solution  which  satisfies  the  appropriate  boundary  conditions  is 
revealed. 

The  detailed  calculation  will  be  published  shortly  in  the 
Journal  of  Geophysical  Research. 
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